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ABSTRACT 

This thesis develops methods of frequency analysis and synthesis 
of digital computer programs describable in the form of a linear difference 
equation with constant coefficients. 

The mainspring of this investigation was the need for dealing 
with control systems consisting of both analog and digital filters. 
Most conventional control systems consist of analog units and operate 
on continuous data, but digital computers use sampled data. A uniform 
treatment of the two types of data is essential in the analysis of control 
systems incorporating a digital computer. !nie conventional method of 
treating systems operating on only continuous data uses Fourier or Laplace 
transformation^ that is, transformation to the frequency domain* The 
conventional method of treating digital programs is numerical analysis, which 
deals almost exclusively in the domain of the independent variable^ that 
is, the time domain. By eaq^loiting and further developing those areas of 
numerical analysis to which frequency-transformation techniques were 
applied, the thesis points the way to a common language of dealing with a 
mixed-data system. 

If data-are sanpled. at equal intervals of time (a practical 
feature), description of a linear computer program always reduces to a 
difference equation. It is possible to describe such a program by a transfer 
function in the frequency domain in a manner analogous to the conventional 
description of analog filters. Whereas components using continuous data 
have transfer functions which are rational functions of the complex frequency 

iii 
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-sT 
variable £, those of a digital program are rational functions of !^ « e , 

where e is the Naperian base and T is the constant interval of sampling* 

Having described the digital coirgjuter with its program by a 
transfer fiinction, one may apply all the techniques of complex-variable 
and transform theory to deal with digital filters. Theorems on realizability> 
stability and other properties of programs are developed, and the amplitude, 
phase and locus of a program are defined. The adaptation of the methods 
of analog filters to digital ones is direct, although the necessary 
modifications are often significant • 

The sjmthesis of computer programs can be conducted along lines 
employed in the synthesis of networks* First, the desired frequency charac- 
teristics of the program are stated^ next, a rational function of 

-sT 
z = e is found which approximates the desired characteristics for real 

frequencies, s » JctJ) finally, the program is realized on basis of the 
approximating transfer function. For facilitating the approximation 
basic entities or blocks of programs are analysed and methods are shown 
by which such programning units can be combined to obtain the frequency 
characteristics of the complete program. Various methods of program 
realization, that is, programming, are developed and compared on the 
basis of time and storage requirements, and criteria are developed to 
permit the choice of the optimum programming procedure by considering 
the mere form of the program transfer function* 

Numerous examples of program analysis and synthesis are shown, 
and one example of synthesizing a program for the compensation of a control 
system is worked out. The latter example shows that the frequency analysis 
of a complete hybrid system can be undertaken along the conventional lines 
and that digital compensation of a control system is possible. 

iv 
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The application of the methods of the thesis to various problems 
in numerical analysis is also shown* The problems of convergence (stability) 
and of truncation errors (approximation) can be analyzed in the frequency 
domain effectively • The study of convergence W conformal mapping is related 
to the usual methods, and a novel way of estimating truncation error is 
shown provided only that the function to which the numerical process is 
applied can be described by its frequency spectrum* 
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INTRODUCTION 

The use of digital computers in control systems is now coming 
into the fore. Unlike most conventional control systems involving analog 
units which operate on continuous data, a control system employing a 
digital computer of the present-day type must use sampled data in the 
part of the system involving the digital computer • Hence, some parts 
of this system use continuous data and others, sampled data. The Fourier 
and Laplace transform methods of' analyzing continuous-data control systems 
is well-known and developed, but the conventional treatment of digital 
computer programs is by numerical analysis or in the time domain* There- 
fore, in order to apply the methods of frequency analysis to control 
systems involving digital computers (mixed-data systems), the sampled 
data part of the system must be described in tiie frequency domain. Some 
work along these lines has been done but it must be furliier developed. 

An analog system is a physical model of a set of differential 
equations; whereas, a digital system is a physical model of a set of 
difference equations. Operational and transform methods have been applied 
to difference equations for some time. In 19hZ Gardner and Banies 
presented a comprehensive and systematic treatment of the solution of 
linear difference equations with constant coefficients by the Laplace 
transform method. However, they do not deal with stability and errors 



^ 

Gardner and Barnes, Transients in Linear Systems , John Wiley and Sons 
New York, 19U2, Chapter IX. 
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which are important in control applications. The control point of view 

is stressed in Tustin's work on time sequences. In 19k9 and 19^0 TustLn's 

2 
method was further developed by Madwed , who shows the relations of his 

aspects of stability, but they do not analyze the errors associated with 

their approximations, 

3 
In the meantime, Hurewicz pioneered the analysis of pulsed 

filters in the frequency domain, developed stability criteria, and 
showed several examples of choosing parameters. It should be noted, 
however, that Hurewicz »s filters are only simple units such as differ- 
entiators and lead networks, which are incapable of performing involved 
computations as a computer can. Also, Hurewicz evaluates the output of a 
pulsed filter at the sampling instants only. The behavior of the filter 
between pulses remains a separate problem, and no reacfy method is pre- 
sented to investigate the whole question in the frequency domain, 

k 
W, K, Linvill shows that sampling a continuous function is. 

equivalent to the modulation of a series of unit impulses by the function. 

The result is a new time function which can be thought of as being applied 

to the sampled data part of the system. Furthermore, this new time 

function has a Laplace transform! thus a frequency-domain analysis is 

possible, Linvill shows that reconversion from discontinuous to continuous 



Tustin, A Method of Analysing the Behavior of Lj.near Systems in 
Terms of Time Series , J.I,£E. Vol, 9it, Part 2A,#1, pp, 130 - lit2, 

2 

Madwed, Numlper Series Method of Solving Linear and Mon-Linear 

Differential Equations , SC,P» thesis in Mechanical Engineering, MIT. 

3 

Hurewicz, Filters and Servo Systems with Pulsed Data , Chapter 5 of 

James, Nichols and Phillips, Theory of Servomechanisms , 

Linvill, ¥.K,, Analysis and Design of Sampled-Data Control Systems , 
Digital Computer Laboratory, MIT, Report R=170, 
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data is a filtering process and also shows what happens when the loop 
is closed on a mixed data system. He is concerned only with the effect 
of sampling on the system and does not consider the influence of digital 
computer operations on the system. 

This report is a sumimry of the work done fcy Salzer • His 
results permit the analysis of linear digital computer programs in the 
frequency domain j i.e., the operation of a digital computer program is 
described by a transfer functiouo Thus the field is opened for the 
complete analysis and synthesis, wholly in the frequency domain, of control 
systems employing digital computers. 

From the frequency-domain point of view, conditions governing 
the realizability of program transfer functions are developed, the problem 
of stability is stadied, and conditions to insure stability are given. 
Three methods of realization of programs from their transfer functions are 
presented, and the time and storage requirements of each are studied. An 
elementary example of transfer function synthesis is given. As in the case 
of network theory, the analysis of a computer program in iiie frequency domain 
is straightforward with a unique result, but the synthesis of a transfer 
function has many alternate realizations. Also as in network theory, the 
characteristics of the transfer function to be realized may not be given 
directly in a foim leading to immediate realization but an intermediate 
approximation problem may need to be solved. The background for solving 
the approximation problem has been set up in that conditions of physical 



Salzer, J.M., Treatment of Digital Control Systems and Numerical Processes 
in the Frequency Domain, SG.D. thesis in Electrical Engineering 
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realizability have been derived and methods of realization of all 
realizable transfer functions have Taretem obtained. While some work has 
been done directly on the approximation problem, much remains to be done 
in this respect. 

The function of this report is to provide a concise picture 
of the frequency analysis of digital control systems and numerical pro- 
cesses • The first chapter describes the processes of sampling and de- 
sampling continuous functions and indicates that sampling is analogous 
to impulse modulation while desarnpling is analagous to ripple filtering 
in demodulation, lliinking of sampling as impulse modulation allows one 
to relate the san5)led to the continuous function in either the frequency 
domain or the time domain. Furthermore, thinking of sampled functions 
as impulse modulated functions allows one to characterize linear computer 
operations on the sampled functions by transfer functions. 

Chapter II derives the conditions of physical realizability 
for computer-program transfer functions, discusses stability conditions 
on these transfer functions, and presents procedures for plotting transfer 

JLOCX # 

Chapter III deals with techniques for realization of transfer 
f\mctions with some attention to the approximation problem, while Chapter 
IV deals with frequency analysis of some numerical integration formulas. 
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CHAPTER I 
DESCRIPTION OF THE SiiMPUNG PROCESS 

1.1 Ssugipling a ContinuoTis Function 

A digital computer operates on numbers that represent samples of 
continuous signals taken at discrete instants of time. The time interval, 
T, between samples is a constant as showi in Figure 1.1, page 7 • In this 
case, the input to the computer is the sampled function, i (t). The com- 
puter senses the amplitude of each of these pulses (as a number) and 
operates on the number. 

The purpose of this chapter is to describe the sampling process, 
to characterize it mathematically, to evaluate how well a continuous signal 
may be represented by its samples, and to show how and under what conditions 
a continuous signal may be recovered from its samples. 

The mathematical model of the sampling process which will be de- 
rived later is very similar to actual physical processes. For exairple, 
assume that i (t) is the voltage across a pair of terminals of some net- 
work. How might it be sampled? The voltage may be sampled by connecting 
a condenser across the terminals, allowing a current flow to build up a 
charge on the condenser until the condenser voltage is equal to the terminal 
voltage, and then disconnecting the condenser. In order that the condenser 
voltage be equal to the terminal voltage at some instant of time, the 
sampling time should be as small as possible. It may be made very small, but 
not zero. The total charge on the condenser is the integral of the curr^t 
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flowing over the time required to take the sample. Thus, as the 
sampling time decreases, the current intensity must increase. 
Physically this is how sampling might be done. Ideally, however, 
we wish to take the sample instantaneously or in zero time. Therefore, 
for ideal sampling in the above example the current flow must be infinite 
for zero time at each sampling instant. Thus, in the ideal case the 
charging ciirrent is an impulse whose area equals the amount of charge 
required to build up the condenser to the sampled -value. Physically, 
ideal sampling is not possible, but the idea permits us to set up a 
model of sampling that can be treated mathematically. 

1»2 Equivalent Mathematical Model of Ideal Sampling - Impulse Modtilation 

The ideal situation in the above example is to transfer to the 
plates of the condenser a portion of charge in zero time, or to "hit" the 
condenser with an impulse of current. The same end can be obtained if we 
modulate the voltage waveform with an infinite series of unit impiilses 
separated by equal intervals, T, as shown in Figure 2. Ttie area of amy 
one of the modulated impulses equals the value of the input fimction at the 

corresponding instant of time. Thus, impulse modulation is ana30.gous to 

1 

the process of sampling. The samples of Figure 1.1 have finite height, zero 

width, and zero area; therefore, the sampled function does not have a Laplace 
transform. The impulses of Figure 1.2 have infinite height, zero width. 



i ^ ' 

The bar (— ) over ift) indicates the sampled functions, 

2 

The circumflex (^ over i(t) indicates the impulse - modulated function. 
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but non-zero area; therefore, the impulse -modulated function does have a 
Laplace transform, which is why this mathematical model has been set up, 

1,3 Use of Impulse Modulated Functions in the Analysis of Linear Digital 
Computer Programs 

A digital computer operates on numbers that occur at discrete 
instants of time, i.e. it operates on samples of a continuous function. 
In the previous section it was i^own that for the ideal case, sampling is 
equivalent to impiilse modulation. If we think of the computer as "sensing" 
the amplitude of samples, we may just as easily think of it as "sensing" 
the area of impulses • With this extension or mathematical model, we may 
analyze computer programs by describing the input to the computer as 
impulses instead of samples. Since a sample does not have a Laplace trans- 
form, while an impulse does, the advantage of this extension is immediately 
obvious. In this mathematical model, both input and output are treated 
as impulses, and both have Laplace transforms. In conventional (continuous- 
data) systems, the transfer function is the ratio of the transform of the 
output to that of the input. Since both input and output of computer 
programs (^en treated as impulse-modulated functions) have transforms, 
we may define the transfer function of a linear computer program as the 
ratio of the transform of its output to the transform of its input. In order 
to carry out this analysis, we must have a knowledge of some of the properties 
of impulse -modulated functions, or iupulsed functions. jQie remainder 
of this chapter is devoted to a discussion of some of these more useful 
properties. 
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l^k Laplace Transforms of Impulse-^Modulated Functions 

Oiip analysis of computer programs is restricted to the cases in 
which the time interval between samples is a constant, T, Thus, the 
impulsed fimctLon can be expressed as the product of a continuous input 
function and an infinite string of unit impulses, the interval between 
impulses being T. 

As the following derivation will show, the process of impulse 
modulation may be readily described in the frequency domain. Essentially, 
since the string of unit impulses (which is the carrier) has all harmonics 
of equal amplitude, the impulse modulated wave has an infinite number of 
side-bands rather than Just the two which are present for a sinusoidal 
jjarrier. The method of the derivation is to'-makel ai Fourier- aa^^ of 

the carrier and to associate each side-band of the impulse modulated wave 
with a Fourier component of Ihe carrier. Let i(t) be iiie continuous input 
function and A:z — Zo^^^^ " ^-^^ ^® ^® infinite string of unit impulses. 
jji (x) = unit impulse occurring at x = O.J Then the impulse -modulated input 
function is, . 



i(t) = i(t) ^^ ,^ p.(t - kT). (1-1) 

To find the Laplace transform of (1-1) let us first find the 
complex Fourier series of the string of unit impulses. 

^a_ „ '-=^:^c e^™-^* (1.2) 



k = "'^ 



tiCt - kT) ^ ^ 



A more complete derivation and discussion of the transforms of impulse 
modulated functions is givsn in Eeference 2, page 3. 
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-lo- 



in (1-2), -H^ = 



2 IT ' 

— =5 — . The c _ s are the complex Fourier coefficients. 



Solving for c in the usual manner we have, 
m ' 



ra 



T/2 ^ 

i j I ^^ t^C^-k« e-^'"^-^* dt. (1-3) 



-T/2 



lE^ writing out a few terms of the series, (1-3) becomes,, 
T/2 



m 



m I I ••••••••• 

-t/2 



+ [i(t - T) + (xCt) + \i(t + T) X.,.. 



-jm-Ht 



dt 



Cl-i*) 



Within the range of the integral, the only term inside the bracket of the 

integrand that is non-zero is the term, p.Ct). Thus (1-ii) becomes, 
T/2 

c, = i J ^(t) e-^--^^ dt. ^^^^ 

-T/2 



Because of the unit impulse in the integrand, the value of the integral is 
just e"*' evaluated at t = 0, which is unity. Therefore, 



_ 1 
m T 



(1-6) 



and the Fourier series of a string of unit impulses is. 



T^^ tt(t-kT) = I 



jmJ\it 



(1-7) 



k = -<»^ 



m = -«^ 
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Then the impiilsed fTmction becomes. 



1 it; = -Y- X e (1.8) 



la = - *»«? 



Now take the Laplace transform of the above equation, 
-iCs) = 1 [Tct) ] = L [^ 2lL -^'"'^* J (1-9) 



The indicated summation can be done after the transformation is made. 



i(s) = i y^ L 1 i(t) eJ""-^* 

m = - <*o 



I' 



(1-10) 



A fimdamental theorem in Laplace transform theory leads directly to the 
following resiilt: 

^ (s) =~ y^ I (s+ jmJX.) (1-11) 

m = - -=*3 

Ihus we see that the Laplace transform of an impulsed function is periodic 
having a repetition interval of j.-ri-. 

An important fact about I (s) should be observed from (1-11), It 
is that there is a unique correspondence between I (s) and I (s) if and 
only if the frequency spectrum of iCt), the continuous time function^ lies 
in the range, ^'j'K.^^'%* 1^ "^^ spectrum of iCt) lies outside this 
range, l(s) will specify the spectrum (in the range -~Q^^co<-^) of a 
continuous time ftinction, but this time fimction will differ from the 
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A. Specturm of i(t) 




♦-^ 



B» Spectrum of i(t) 



Figure 1«3 Unique Correspondence Between l(s) and l(s) 
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original time function. Thus, there is a limitation of bandwidth 
caused by sampling. Figure 1,3 illustrates the case of unique corres- 
pondence, and Figure l.ii, the case in -tdiich the spectrum of i(t) is too 
Tijide* 

As given lay (1-11), l(s) consists of an infinite number of termsj 

however, an infinite series is difficult to handle, and it is desirable 

to have a closed form expression for l(s). This can be obtained from 

K. 

the partial fraction expansion of I(s)* Consider a typical term, — , 

s •" s « 
1 

of the partial fraction expansion of I(s)» Referrihg to (1-11) we see 
that corresponding to this typical terBi,^(s) will have a typical series of 
terms of the form. 



h ^^Z^ 1 

T > s - s. + jk J:\^ • 
k= -«» ^ 

Thus we see that the pole at s = s. is repeated -eji infinite number of times 
at intervals of j-^, the line through these poles being parallel to the 
imaginary axis in the s - plane. 

A closed form equivalent of the above typical series can be 
obtained by a change of variable in the following equation. 



n r cot ,T a = 1 + 2 Z- 2 J 1 ^^_^^ 

n = 1 



1 

Knopp, "Theory and Application of Infinite Series", New York, 19U8, p. 1419 
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B, Spectrum of i,(t). 
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C« Spectrum of i^Ct) that would produce the same sairpled ftmction 
as (B)<. 



Figure iJi. llltistration of bandwidth Limitation Caused t^y Sampling< 
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-1^- 



CO 



Divide each side of (1-12) by z", and make the change of variable, ^ = j — . 



Tf cot j 



ff CO __ ,XL, 



Go 



- j2 ^ ^ 



n = 1 



31?'*" ^ 



(1-13) 



Multiply both sidss of (1-13) by j/fu and obtain. 



i JL. cot i ILS = JL coth HL^ = 1 + 



n = 1 



2co 



CO 



2^2' 



+ ttSl, 



(1-lU) 



The infinite series of poles of I(s) corresponding to a pole of 
l(s) at s. can be put into a form that is identical to the right-hand member 
of (l-li^) as follows! separate the terra for k = 0. 



Qear 



k = -^ s - s^+ jk.Jl. 



s - s. 

1 



*1 



so. 



k = 1 



1 1 

s-s. + jk-O- s-s. - jkXU 



^1 



(1-15) 



Combine Ihe two terms in the s-ummation. 



^ 11 

s-sT + JkJL ^ " \ y 



2(s - s.^ ) 



^ 



T^-r 



«»<=» 



k= 1 



(s - s.)^4. k^ iV^ 



SL-16) ) 



A comparison of (l-lii) and (1-16) shows that. 



k = - o» s - s . + j k-ni« 



1^ eoth ^ ^" - ^1 > 



-O- 



_ru 



(1-17) 



Thus we have ihe following closed form equivalent of the typical series of 
I(s), 



K. 



Cko 



J 



k= - oo 



s - s. 4J ^kJ - L = 4 °°^ I ^^ - ^i>> 



(1-ia) 



Report H-225 



-16- 



for a pole of l(s) at s = s. , Therefore, corresponding to the partial 



/Va 



fraction expansions of l(s), we have the following series for l(s). 

n 
iCs) = I \ K^ coth I (s - s^) , 

i = 1 

where "n" is the total number of poles of l(s), taking into accoimt multiple 
poles» 

Let us now investigate the limitations on the positions of the 
poles of Ifs) due to sampling. Consider an infinite strip of width -H^ in 
the s. -plane and parallel to the real axis as shown in Figure 1.5. Assume 
that all the poles of Ifs) lie within this strip and in the left half plane 
(LHP). Thus, iCs) has 



jw 






% s. 



s-plane 
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s*. 
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L^ 
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. -Tu 



Figure 1.5 Infinite Strip Containing Poles of I(s) 
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poles at these points plus poles at points shifted from the s. s and 

s*. 's C* means conjmgate) by the distance « jk XL • Since I(s) has 

poles only in the strips "being considered, there is a one-to-one correspondence 

between the poles of ifs) and l(s) that lie in the same strip* However, 

if l(s) had poles outside this strip, there no longer would be this one-to-one 

correspondence • 
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CHAPTER II 
TRANSFER FUNCTION OF COMPUTER PROGRAMS - REALIZABILITY AND STABILITY 

Using the properties of impulse -modulated functions given in 
Chapter I, we are now ready to investigate transfer functions of computer 
programs. Our interest in program transfer functions is much more than 
academic. The transfer function describes the program completely and 
with it we can analyze and synthesize control systems enploying digital 
computers by conventional frequency domain methods* 

In this chapter a linear digital computer program is defined in 
terms of the mathematical model of sampling set up in Chapter I, its transfer 
function is derived, and methods for deteimining the realizability and 
stability of transfer functions are given. Several examples of stability 
determination are also presented. 
2.1 Transfer Ftinction of Linear, Real-Time, Digital Computer Program 

As pointed out in Chapter I, the input to a digital computer 
may be assumed to be an impulsed function, for purposes of mathematical 
analysis. A linear program of a digital computer operating in real time 
is one in which the present output is a linear function of the present 
and past inputs and the past outpiits. The general form of this relation 

m ^ n 

"oCt) = 5^ aj^ itt - kT) « y^ h^ oCt - kT), (2-1) 

k « k « 1 

in which all a, *s and b, 's are real, and T is the time between samples* 
The time required for the cou5>utation must be less than T if each 
calculation is to be con^leted before the next input arrives. 
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Taking the Laplace transform of (2-1) yields, 

-0(3) .'l(s) y^ a^ e-^^^ - g(s) y^ h^ e"^^^. (2-2) 
k « k = 1 

As in continuous data systems, we will define the transfer function 
of a computer program as the ratio of the transform of the output 
to that of the input. Let WCs) be the transfer function of a computer 
program^ then, 

W(s) =- Bil . (2-3) 

Ks) 

Solving for tJts)/ I(s) from (2-2) we obtflpLn, 





lis) 


■fz -. •-' 


w(s) - 


k « 





{2-W 



as the transfer fxinction of a linear, real-time, digital computer program. 
With the understanding that b = 1, (Sf-ij.) becomes, 

m 

"S -ksT 
-^ a e 

¥(s) » k :» ^ . (2-5) 

^ ksT 



k « 

The inverse steps from (2-5) to (2-1) are uniquei therefore, 
(2-5) is the general form of the transfer function of a realizable, linear, 
digital computer program^ . Thus, to be realizable, the transfer function 

of a linear, digital computer prograa must be expressible as the ratio of 

-sT 
two poljmoraials in e • The criteria for stability will be discussed in 

a later section. 
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It has already been sho-wn that the Laplace transform, Ivs), of 
the impulsed input function is periodic of period •-^^^, as seen in (l-ll)« 
By showing that >r(s) is also periodic with the same period, we can prove 
that 0(s) is also periodic of period Jl^ e A typical terra of either numerator 
or denominator of W(s} contains e o For s — >s + jra-H. (m is a positive or 



negative integer), the typical terra becomes 



-k(s + jmn)T _ -ksT -jkmi^T 

e / — e e « 

As T-Q. = 2ir and k and m are integers^ the second factor is, 

- jkm •CHS -j2irtcra ., 
e s! e «* lo 

Wor^.o «*kfs * JBijfDT ^ -ksT 
Hence, e = e « 

Therefore, the terms of the numerator and denominator |i)^ W(s) are periodic 
of period ,CL^ and so is W(s)« In equation form this means, Wfs) « 
W(s ♦ Jm_Q}, for m a positive or negative integer • The product of two 
periodic functions is also periodic o Since 0(s) « W(s) l(s), llTCs) i? also 
periodic of period./!., as indeed it should because the computer output is 
also sampled. 

Since all the coefficients of (2«^) are real, it is readily 
seen that W(s) « ¥(s ) , in which the asterisk means conjugate* For real 
frequencies this becomes ¥( joa) = W(=jo5) « This fact together with the 
periodicity of W(s) tells us that W(s) is completely specified for all s 
if it is defined over the range, ^m^-^o 

Stimmarys In order to be realizable^, the transfer function of 
a linear digital coi^uter program must be expressible as the ratio of 
polynomials in e° «. W{s) is periodic of period^TLl ioeo^ W(s) « 
W(s ♦ jm-/l)« Specification of W(s) over the range, O^os^-fe completely 
determines WCs)o 
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2<»2 Stability of Programs 

¥e have expressed the transfer function of computer programs as 
a function of the complex frequency **s**| therefore, the same methods of 
investigating stability as used in network analysis and servomechanisms 
are applicable <. The general necessary and sufficient criterion for stability 
of a \init is that its transfer function have no poles in the right half 
s-plane (HHP) or multiple poles on the jos^axiso In network analysis the 
frequency-'domain method used to study stability is to map a contour 
enclosing the right half of the s=plane (the contour is usually the jco-axis 
and an infinite semicircle) into the W-plane. Because of the transcendental 
nature of the transfer function of a realizable computer program, the 
mapping contour in the S'-plane must be modified** 

As we have shown before, the transfer fiinction of the computer 

program is^ 

^ "^ -ksT 

¥(s) « P(s) « k « (2-6) 



Q(s) ^5— ^ «=ksT 



y b, e 



k « 

in which P(s) is the numerator and Q(s)^ the denominator of W(s)| and it is 
assumed that P(s) and Q(s) have no common factor* Both P(s) and Q(s) are 
entire transcendental functions having as their only singularity an essential 
singularity at infinity.,* Hence, we see that the only singularities of 
¥(s) in the finite s-plane are poles, and these poles occur at the zeros 
of Q(s)o Our stability criterion is that there be no poles of W(s) in the 
RHP and only sipple poles on the imaginary axiso Therefore, in order for 



* For a further discussion of entire transcendental fxinctions, consult Knopp, 
^Theory of Functions," or Guillemin, **The Mathematics of Circuit Analysis*** 



Report R-22^ <»22-> 

the program to be stable ^ QCs) must have no zeros in the RHP and only simple 
zeros on the jco«»axiso To investigate the possibility of Q(s) having zeros 
in the RHP or on the imaginary axis, we may take advantage of the periodicity 
of Q(s)o In proving that "W"(s) is periodic, it was shown that e~ is 
periodic property* Therefore^ if Q(s) has a zero in the RHP, it must have 
one in the semi«inifinite strip shown in Fig* 2 do 



s -plane 



^c- 




cr ->^ 



sT 



Figure 2«1 Semi°inifinite strip of s~plane that must have a zero 
of Q(s) if Q(s) has any zeros in the RHP, 

"-sT 
Consider the map of the contour of Figo 2ol into the e plane 

Let us begin the path at the origin in the s=plane and encircle the strip 

in a clockwise direction, corresponding to increasing frequency^ It is 

readily iinderstood that corresponding path and enclosed region in the e 

plane is shown in Fig» 2»2o The origin of the s~plane maps into the 

point (1,0) in the e" plane « The corresponding sections of the path are 

marked by small letters on both contours » In Fig<, 2o2 we see that the 

paths (b) and (d) cancel leaving the annular ring as the region confomal 

to the strip of the s-plane that is under consideration. As O; (of 

Figo 2«l} approaches 00, the radius of the circular path (c) in Fig* 2,2 

approaches zero. Thus, the conformal map of the indicated strip consists 

of two separate contours? one, a unit circle centered at the origin 
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-sT 
e -plane 



unit circle 



Figure 2.2 Conformal map of the semi-inifinite strip of Fig» 2*1 



into the e 



'St 



plane. 



and the other an infinite simally small circle that excludes the origin 

in this particular case» Only a slight extension of the foregoing procedure 

°=sT ""ksT 

is required to deteraiine the map of powers of e • The map of e 

-sT 
will appear like that of e except that each of the two separate paths 

will be traversed "k** times| the region excluded by the infinite simally 

small circle will be that at the origin o Thus we see that the map of 

this serai -infinite strip of the s«plane is effecti?ve in handling the 

essential singularity of Q(s) at o© « 

Now, consider the conformal map^of the serai-infinite strip of 

Fig. 2»1 into the Q-plane, Remembering that Q(s3 « Ifb, e" **" ^o® * 

■»nsT 
.•• ■*• b e , we see that the map of this strip into the Q-plane will 
n 

exclude the point (1,0), (the map of each term except the first excludes 
the origin). This eliminates the need for mapping path (c)o Moreover, 
since the paths (b) and (d) cancel, we meed to plot only the paths (a) 
and (e)« In other words the only part of the s-plane contour that we 
need to plot in order to deteimine the locus of Q(s) is the part of the 
contour that lies on the imaginary axis. This contour in the Q-plane 
Wi31 encircle the origin z-Nitimesi^ the counterclockwise directi^n^ where 
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z is the number of zeros and N is the number of poles of Q(s) (taking into 
account their multiplicity) in this strip of the s-plane* It has already 
been pointed out that Q(s) is an entire transcendental function and, therefore, 
has no poles in this strip. So, N = 0, and the contour in the Q-plane will 
encircle the origin Z times (clockwise is to be understood). The condition 
for stability of W(s) states that Q(s) must not have any zeros in the RHP 
or any multiple order zeros on the imaginary axis. Therefore, the map in 
the Q-'plane must not enclose the origin ! Z must be zero. If Q(s) has zeros 
on the imaginary axis, the Q-plane locus will pass through the origin. In 
this case, we must determine the order of the zero. The following method 
can be useds Asstime that the locus in the Q-plane passes through the origin 
for s * jco. * Then Q(s) must contain the factor (e~ - e"*^ i ) ^ere 
n is the order of the zero. Divide Q(s) by (e" - e""*^ i ) • If there is 
no remainder, the zero is of higher order than the first and the program 
will be unstable. 

In addition to determining the stability of programs, conforraal 
maps give an indication of the degree of stability or instability and an 
approximate value of the frequency at which the program is or may become 
unstable* The amount by which the locus in the Q-plane misses encircling 
the origin gives a measure of the stability of the program. The farther the 
locus is from the origin, the more stable or convergent the program. The 
frequency corresponding to the point on the Q-plane locus nearest the origin 
is approximately the frequency at which the program is or may become unstable, 
or at which it will oscillate in a damped fashion. 
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In addition to expressing the program transfer fimction as a function 

-sT 
of "s** we may also write it as a fiinction of e « Make the change of variable, 

z = e°° o Then we may define a new fimction, 
VCz) 



NCz) . 


> V 

k - 


D(z) 


> ..V 



(2-7) 



k = 

/ 2 Ifl 

rrt \ SI + a- 2 ♦ a«z ♦ ««<ioo<>oo« + a z 
V(z} - o 1 2 m 

1 + b. z + b-z ♦ ,..•.««,. + b z 
12 n 



C2-8) 



It is readily seen that the right half of the s-plane maps into the inside 
of a unit circle centered at the origin in the z-plane» The imaginary axis 
of the s^plane becomes the unit circle in the z-plane (see Fig, 2»3)« 

s=plauie 





Figure 2o3 Map of right half of s -plane into z-plane 
Therefore^ if the program is to be stable, all the zeros '6f D(z) must lie 
outside the unit circle except that single order zeros may occur on the 
unit circle o In other words, the magnitude of the roots of D(z) must be 
greater than or equal to tanity, and the roots of unity magnitude must be 
simple* 
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Sxunmary: To test for the stability of a prograia, map the semi -infinite 
strip of Fig. 2.1 into the Q-plane £q(s) U the depgrninator of the program 
transfer function^ If the locus in the Q-plane does not enclose the origin, 
the program is stable or convergent. If the locus passes through the origin, 
Q(s) has a zero and the order of this zero must be determined* If the zero 

is of first order, the program is stable j otherwise, unstable. An alternate 

-sT 
method is 8 Make the change of variable, z « e , and find the magnitude 

of the z-roots. If each root has either a magnitude greater than unity or 

equal to unity and is simple, the program is stable or convergent^ othemise 

unstable or divergent. 

2.3 Loci of Q(s) 

In the previous section it was demonstrated that the stability 
of a program can be determined by mapping the contour enclosing the 
semi -infinite strip of Fig. 2.1 into the Q-plane. It was also shown that 
the only part of this contour that we need to plot is that on the imaginary 
axis. The paths (b) and (d) cancel and the path (c) excludes the point 
(l,0) in the Q-plane,. Hence, we are interested in the properties of Q(jcd) 
and its locus in the range, - J'^&i ~?' 

QCjco) has several properties that are helpful in determining its 
locus • 

(1) QCjGj) is periodic of period-^. This was proved in the 
previous section. 

(2) Q( j(a) » Q(-;3co) . This property follows directly from 
the fact that Q(joo) is a polynomial in e""*^ , and as a 
consequence, the locus of QCjco) must be symmetrical about 
the real axis. 
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(3) At ca « and ♦* Qfjco) is real and its locus at these two 

points crosses the real axis either normally or tangentially* 
This statement is proved and elaborated upon in Appendix A* 

As a consequence of the first two properties, the locus of QCjco) 
for ^ CO ^'Q' completely determines the locus in the Q-plane • The locus 
for -"x*:^© ^0 is just the mirror of that for the positive values of o). 
Thus the first two properties result in a substantial reduction in the amount 
of work required to plot the locus of Q(^(»). The third property enables 
one to determine accurately the shape of the locus in the neighborhood 
of 00 « o and ♦ ^ * 

Several methods may be used to determine the locus of Q(jco)» 
Three of these are: (l) add the loci of the individual terms of the 
polynomial (each locus is a circle) to obtain that of Q(ja>)^ (2) factor 
Q(jco) and multiply the loci of the factors; and (3) express Q(jco) in the 
form R(ca) /0(<a) and make a point by point plot. The method that is best 
to use depends on the particular Q(jco)« However, it is to be expected that 
the extra analytic work required in methods (2) and (3) will result in 
less graphical work and more accurate loci* .None of thp methods will be 
discussed, but t^ey will be illustrated* 

Let us now consider several exarnples of loci of Q(Jcd)« 

1, Let Q(s) « 1 - 0.8 e"®*^ ♦ 0,3 e""^^^ (2»9) 

Take the derivative with respect to s« 



ds 

^ « Oo2Te"^^ {k - 3e"^^) (b) 

OS 



(a-10) 
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For neither s « nor — j4 is the derivative zero, so the locus is normal 
to the real axis at both points* Fig» 2»k (a) shows the locus of Qfjco) and 
how it was obtained (for the point cot « r— ) from the loci of the individual 
terms • The locus for negative values of co is shown by the dashed curve, 
since it may be obtained from the other half of the locus. After this, only 
the locus for positive co will be drawn* The locus does not enclose the origin^ 
therefore, a program whose transfer function has the denominator, 
1 - O.Se*"®^ ♦ 0,3e"^®^, will be stable* 

2. Let Q(s) « 1 - Oo8e"^^ ♦ OoUe"^®^ (2-11) 

Take the derivative with respect to s« 

~|| = o.BTe"^'^ - 0.8Te"^^^ (a) 

- O.STe'"^'^ (1 - e""^^) , (b) 



(2-12) 



For s « 0, the derivative is zero^ Q(s) has a saddle point here* Q(s) can 
be rewritten as, 

Q(s) - 0*6 * 0*i4 (1 - e*"^*^) (2-13) 

which brings the saddle point into -evidence* In this case p = 2, so at 
CO « 0. the locus is tangent to the real axis. At s « ♦j^j -r- / 0, so the 
locus is normal to the real axis at this" point. Fig. 2*1^. (b) shows the 
resulting locus. It does not enclose or pass through the origin^ therefore, 
this QCs) will lead to a stable program* 

3* Let Q(s) - 1 - 0*8e'®'^ + O.^e"^^'^ (2-lii) 

Take the derivative with respect to s« 

IS - o.8Te-=^ . Te-2«^ (a) 

- 0.2Te-^^ ik - Se'^n (b) 



Report R-225 ■°29<=> 

For neither s « nor *i^±s the derivative zero, so the locus of Q( jco) is 
perpendicular to the real axis at both points. The locus (as shown in 
Fig. 2,k (c)) does not enclose the orig:|.nj therefore, this Q(s) is the 
denominator of a stable program transfer function* 




(a) Q(s) « 1^.8e"*^ ♦ Q^"^* 




(b) (}(s) « l-0.8e'^^ ♦ OJie"^^^ 



M Q(s) - l-0,8e-^^ ♦^e"^^^ 



Q-plane 




(P 
o 

PO 



o 
t 



Figure 2.k Loci of some typical Q(s) 
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CHiLPTER III 

ANALYSIS AND SYNTHESIS OF LINEAR, DIGITAL COMPUTER PROGRAMS IN THE 

FREQUENCY DOMAIN 

In the first part of this chapter the analysis of transfer ftinctions 
is dealt with by expanding the transfer function into partial fractions. 
Next, programs are realized from transfer fiinctions by three methods: direct 
programming, cascade programming, and parallel, prograraraingj and the storage 
and time requirements of each are presented* In the last part of the chapter 
a short, general discussion of synthesis is given, and one possible synthesis 
procedure is illustrated by the sjmthesis of a program for differentiation. 
3»1 Response of Programs at Real Frequencies 

The input to a computer has a certain frequency spectnam, and 
in order to analyze the action of a computer program on this input function, 
we need to have a knowledge of the frequency response of the program. Thus, 
we are interested in the locus of ¥(jco), the map of the jco-axis of the s-plane 
into the W-plane. A familiarity with the frequency characteristics of the 
simple transfer functions is essential for the understanding of the possibilities 
and limitations of more complicated ones. 

In many cases the desired locus of the transfer function of a 
digital computer program is given, and the problem is to approximate this 

locus by that of a realizable program^ i.e., by a ratio of polynomials in 

— sT J- \ 

e • A study of the loci of typical terms of ¥(;5q)) is helpful in making 

this approximation. 

3*2 Analysis of Building Blocks of Transfer Functions 

As we have seen, the transfer function of a linear, digital 

computer program is most generally expressed as the ratio of polynomials in 
-sT 
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«sT ^ „2sT ^ ^ -msT 

WCs) = 2i£l = S—1 2 m (3^^) 

I(s) ^ ^ , -sT ^ , -2sT ^ ^ «-nsT 
1 + b.e + b^e 4- o,ooo + a e 
12 n 

A partial fraction expansion of ¥(s) can be made, and we may call the individual 
terms of the expansion the basic building blocks of a program transfer functiono 
In general W(s) may be broken up into a polynomial plus first and second degree 
partial fractions (from the real and conjugate complex roots of the denominator, 
respectively) <> 

In analysing computer programs in the frequency domain [finding 
the locus of W( jco) / several methods can be usedo Two of these ares (l) 
Find the loci of the numerator and denominator pol3momials and then divide^ 
(2) expand ¥(s) into partial fractions, find the locus of each term of the 
expansion, and add the resultant locio In most cases the first method 
is easier to use, but the second is included here because of its connection 
to the synthesis of program transfer functions (approximation of a desired 
locus by a sum of the basic building blocks)* A familiarity with some of 
the possible loci of polynomials and first and second degree partial fractions 
is an aid in the synthesis procedure* 

The loci of second degree polynomials have already been discussed, 
and the locus of a fourth degree polynomial will be illustrated in connection 
with polynomial building blocks© Since there is only a short step from 
the loci of polynomials to the locus of a transfer function, the first 
method of finding the locus of W(jco) will not be discussedo 

For use of the second method, we will investigate the loci of typical 
terras of the partial fraction eacpansion of ¥(jco)o 
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3*21 Palynomials 

A typical pol'ynoinial transfer ftmction is of the form, 

W(s) . ZZ V"""'"^ 0.2) 

For s «■ j^co, the locus of each term of the polynomial is a circle • For 
Q(s)> the constant term is unity, but for polynomial building blocks, the 

constant term may have any real value. In the previous tjhapter, three 

-sT 
examples of the locus of second order polynomials iii e are given, so now 

let us find the locus of a fourth order polynomial. Let, 

W(s) - i^ (13 ♦ Ue-=^ - 36-2='^ ♦ Se'^"^ - e-^"') (3-3) 

First examine the function for saddle points. 

^ . I5 UTe-^T ♦ 6Te-2^'f - 6Te-3^% l,Te-^«') (3-U) 

§ ..^ e-^'" (2 - 36-"'' ♦ 3e"2'" . Ze'^'^). (3-5) 

The derivative is zero for s « 0, therefore ¥(s) has a saddle point there • 
W(s) can be written in the form 

W(s) - 1 . 1^ (2 ♦ e-2^^) (1 - e-^'^) , (3-6) 

which brings the saddle point at s « G into evidence. The saddle point 
is of first orderi therefore, the locus is tangent to the real axis at 

s « 0. W(jQj) is 

2 

^(jfio) - 1 - I^ (2 * e"^^^) (1 - e"^"^) (3-7) 

In this case it is easier to determine the locus of W( jco) by 
plotting the second term of (3-7) and shifting the origin one unit to the 
left. A convenient way to find the locus of the second term is to let 

- (2 + e-i^^t (1 - e-^"*^)^ - Re^^ , (3-8) 
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where both R and are functions of coTo Some trigonometric manipulation 
yields, 

R « 2(1 - cos coT) ]/$ * h cos 2coT (3-9) 

and 

(3-10) 



. ^^-i Ptan a)T(5 > tanV)l 
L tan mT - -^ J 



The resulting locus of W(ja)) is shown in Fig, 3«1» 



>< 



¥(s) « \^ ij (2 ♦ e"^^^) (1 - e"^^) 



<■ 



TT 10,2 




4 



* -t- i 1 4 



<aT - 




Figure 3.1 Locus of a Fourth-Degree Polynomial T ransfer Function 

3»22 First-Degree Partial Fractions (real roots) 

-sT 

In the partial fraction expansion of a rational function e , 

a typical terra has the form, 

^ . (3-11) 



W^(s) 



1 4. p e 



»sT 
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If the constants ©^ and ^ are real, then (3-11) cai be considered 
a basic building block. In this section we will consider the case in which 
o<^ and p are real. The case of conplex constaits is considered in the next 
section. 

First, we may set c>C m i because it is merely a scale factor* 
Second, the magnitude of p must not be greater than unity for W-(s) is 
then unstable. If |p[ ^ 1, the typical term is stable. 

To determine the loci of typical terras of the form (3-11 ), we need 
only apply some of the rules of the loci of con^lex functions. First, note 
that the locus of 1 ♦ p e"*' is a circle of radius |p| centei*ed at the 
point (l.O). To find the locus of W^(jo)), the inverse of a circle must be 
found. If )p| "1, the circle (1 + e"*' ) passes through the origin, 
and its inverse is a straight line parallel to the imaginary axis. If 
/p|^l, the circle does not pass through the origin, and its inverse is 
another circle. The loci of two stable first-degree partial fractions are 
shown in Pig. 3*2. 




b» W(s) « 
Figure 3»2 Loci of First-Degree Partial Fractions 



1 •»• e 



-sT 
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3 023 Second-Degree Partial Fractions (complex roots) 
If a typical term, 

^M) « ^ m , (3-12) 

— sT 

in the partial fraction expansion of a rational function of e has complex 



constants, there will be another term, 

-sT 



M (s) . -3 ^ (3.13) 



1 ♦ p*e 

in the expansion whose constants are the conjugates of those of W-(s)» 

(The asterisk means conjugate •) Both W_(s) and W«(s) will be stable if and 

only if Ipl^l. If the rational function has real coefficients (as in 

practical problems) the terms such as ¥.(s) and W^(s) must come in pairs • 

Add the two and obtain a typical second degree partial fraction with real 

coefficients. 

*\ ^ "^ \ —sT 

In this section the discussion is restricted to second-degree partial fractions 

~sT 
whose denominators have complex roots of e « To siR5)lify the analysis, 

¥o(s) can be written as, 

1 ""SX 
+ a.e 

1 + b-,e ♦ D_e 
1 2 

which differs from (3-lU) by only a constant multiplying factor. 

To insure that the roots of the denominator of (3-15') are complex, 

2 
b^ ^ l|.b^» ¥ith this condition imposed on the constants of the denominator 

of (3-15), it can be considered a basic building block of a program transfer 

function. 
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Now let us consider the stability of such a building block* 
¥^(s) will be stable if and only if both ¥. (s) and W^(s) are stable. A 
comparison of (3-lli) and (3-15) shows that, 

b^ - p 4- f, and bg « pp* « |p(^. (3-l6) 

Therefore, the necessary and sufficient condition for the stability of 
WmCs) is that, 0^ bjj ^1* We may combine this with the condition for 
complex roots in the denominator of W-(s) and obtain, 

[t] < ^2<> ^^"^"^^ 

as the necessary and sufficient condition that insures stability of W«(s) 
and complex roots of its denominator. 

In Fig. 3»3 there are plotted three loci of building blocks 
of the form (3-15) • All three are stable. It sho\ild be observed that by 
changing only the numerator of the transfer function, three completely 
different loci have been obtained. 

By adding building blocks of the form discussed in these section^, 
a desired frequency characteristic can be synthesized. The synthesis of 
a differentiating program is discussed in a subsequent section. 




1 - O^Se"^'^ 
1 - 0.8e"^^ + O.Ue"^^^ 




^ 



floT = 



1 - O.Ue 



-sT 



1 - 0.8e'^^ * O.Ue"^^^ 




-sT^ 



1 - 0.8e ^^ ♦ 0*Ue ^^^ 



Figuire 3*3 Loci of Second-Degree Partial Fractions 
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3 #3 Realization of Programs From Their Transfer Functions 

The purpose of this section is to develop and compare methods by 

-sT 

which a program can be derived from a rational function z « e which is 

the transfer function of the program* How this rational function is arrived 
at in the first place is the concern of the last section of this chapter* 
3*31 General Considerations in Program Realization 

III choosing a particular method of programming, one may consider 
the following factors: storage requirements and time requirements. To a 
certain extent one of these requirements can be reduced at the expense 
of increasing the other and the optimum method will depend on the particular 
application* It is necessary, therefore, to make available various possible 
methods of programming and to form some idea about the requirements of eachj 
intelligent program realization can then be adapted to each application* 

In the consideration of storage requirements of linear programs, 
it is convenient to distinguish three tj^pes of storages data storage, constant 
storage, and instruction storage • The data are the successive sampled 
values of input and output. The complexity of a program is closely related 
to the nuinber of constants and to the age of the data to which the program 
refers. The program can be divided into arithmetic and manipulative parts* 
The number of arithmetic operations involved is roughly proportional to 
the number of constants, each implying a multiplcation (of a piece of data 
by the constant) and an addition (of the product to the other terras). 
The number of manipulative operations is related to the "age" of the oldest 



1 

It is understood that in a general-purpose computer there is no physical 
difference between the storage re^sters containing numbers or instructions, 
and any register may hold either kind of information. The distinction 
made here is only for the puipose of discussion* 
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data used, where age is escpressed in terms of sanpling intervals* All 
"yoTinger** data must also be stored even if not used at each calculation 
(the corresponding constants being zero), for eventually they will become 
the oldest data. After each calculation of a new output value, the manipulative 
instructions shift each piece of data to a storage location at which an older 
piece of data has been, the oldest data being lost. The manipulative program 
is seen to rearrange the data storage in such a manner that at the new sampling 
point the same arithmetic program will calculate a new output value • 

The time requirement of a program is the product of the number of 
instructions to be carried out and the average duration of an instruction. 
The latter factor depends on the physical characteristics of a particular 
computer and is more or less fixed| the number of instructions performed 
in sequence, however, depends in part on the manner of program realization. 
In each particular realization a significant trading of time for storage 
is possible by so-called cyclic procedures* One notes that often the calculation 
of each term in a program involves the same sequence of arithmetic operations. 
The siii?)lest and fastest procedure is to store as many of these sequences 
as there are tenns to be calculated o Considerable storage may be saved, 
however, by storing these instructions only once and cycling through them as 
many times as there are terms to calculate <, Unfortunately, the time requirement 
increases considerably, for in each cycle the addresses of the instructions 
must be adjusted to make them refer to different storage locations for different 
terms and the number of cycles must be counted to permit termination of 
the cycling process. 

The following sections and related appendices will serve as specific 
illustrations of these considerations in programming. 
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3*32 Direct Regression or Direct Programming 

The starting point of onr realization procedure is the general 



expression for the transfer function of a linear program, 

-sT ^ «2aT^ ^ -msT 
a + a,e + a^e •»'•.•• ♦ a e 
o 1 ^ m 



W{s) - . g^T .„3. r . (3-18) 

1 ♦ D_ e + Drt© ♦ • • • + D e 
12 n 

In order to interpret a program in the time domain it is necessary to eliminate 
fractional expressions. The most straightforward way of doing this follows 
directly from (3-l8)« From it we can obtain 

^(s) « Ca^ ♦ <.,• -^ a e"^^^) l(s) - (b,e""^^ * ... 4r.b e"^^"^) "Sts). 

u m L • n 

(3-19) 
The inverse transform of this expression is 

^Ct) - a^Ct) ♦ a^Itt » T) ♦ „.. + a^(t - mT) - b£o(t - T) 

" ... ^ h^it - nT), (3-20) 

where 'o(t) andlCt) are inpulse-modulated (sampled) time functions having 

ibhe value zero everywhere except at the sampling points* In terms of some 

continuous functions b(t) and i(t), which agree with the area-values of 

'o(t) and i(t) at the sampling points, (3-20) is often written as 

o.«ai. + a_i.- + «oo + ai. -b^o,--»i*-bo. , (3-2l) 
J o J ^ a-1 m j-m 1, 0-1 n J-n^ 

where j signifies particular sampling point and j-k the k-th preceding sampling 

point. Eq. (3-2l) is more familiar to the numerical analyst than (3-20), 

but the two are entirely equivalent and are called regression formulas. 

These equations state that the present result (output) is computed by a 

finite linear combination of the present and past input values and of 

past results (output values). 

Several characteristics of regression formulas should be observed* 

If the right side of (U-20) or (ii=2l) has at least one non-zero b, , then 
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the present output depends on at least one previous output, which in turn 
depends on an output further back and so on. It follows that the present 
output value is affected by output values as far back as the start of the 
problem and therefore, also by input values that far back. Thus regressing 
to a finite number of output values corresponds to regressing to an unlimited 
number of input values. This aspect of the regression equations is important 
and will be further enq:>hasized in the following sections* 

Interesting conclusions can be drawn concerning memory requirements 
on the digital computer by considering the actual programming of the regression 
formula, (5-21). Assuming that none of the coefficients a, and b are zero, 
one can easily see that in order to calculate a new output value o ., when 
a new input value i. is received, m previous input values and n previous 
output values will have to have been remembered, requiring m + n memory 
positions for data where m and n are the subscripts of the last non-zero 
coefficients, and furthermore, it is necessary to store all these data even 
if some other coefficients are zero because at the next sampling point 
the same pieces of data will be associated with different coefficients. 

It can be stated that, at least when programming is done by the 
illustrated direct regression metho'^j, the data memory consists of m ♦ n 

registers (memory positions) where m and n are the degrees of the numerator 

"sT t \ 

and denominator polynomials in z =. e of the program transfer function ¥(s). 

Actually this data storage requirement may be reduced, as will be shown in 
Sections 3.33 and 3»3U« 

To be able to make comparisons between the various synthesis 
procedures it is necessary to do the actual programming. This exercise 
is left to the appendices, and the results will be con^ared after the 
other synthesis procedures will have been discussed. In Appendix B the 
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arithmetic and manipulative parts of the direct regression program are first 
constructed separately; then a new more compact program is shown which inter- 
leaves the arithmetic and manipulative instnactionso Although the Whirlwind 
code is used, the results and conclusions can be considered quite general 
in view of the fact that the instruction complements of most general-purpose 
digital computers are conspicuously similar* 

3 •33 Cascade Programming 

-sT 
If the mimerator and denominator poljmomials in z « e are 

factored, {U-l8) takes the form 

(1 ♦ c,e-^^) (1 + c^e"^'^) •.. (1 ♦ c e"^^) 

W(s) - a^ ^^—m ^—zffi ^S^-Tsr- , (3-22) 

° (1 ♦ dj^e"^^) (1 + dge"^^) .,. (1 * d^d"^^) 

where -(l/c, ) and -{l/d, ) are the roots of numerator and denominator 
respectively, when considered as polynomials in z» Because the coefficients 
a, and b, of these polynomials are real, the c, and d^ will also be real 
or will come in conjugate pairs* At any rate, it is possible to group the 
monic factors of (3-22) in some manner 

¥(s) = ¥^(s) ¥^(3) •.. Wp(s), (3-23) 

where each W^(s) is of a rational form in z having a numerator and a 
denominator of not higher degree in z than ¥(s) itself has. 

The form of (3-23) reminds one of the transfer function of cascaded 
linear units in a servo systemo Cascading means that the output of one 
unit becomes the input to the next one* There is no difficulty in using 
the same interpretation to define cascaded programs o At every sanpling 
point the output of each regression equation is used in calculating the 
output of the next one* 
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To be more specific, let us assiirae that? (l) W(s) is a proper 
rational fraction in z, that is, m < n ^ (2) all roots ~l/c^ and -l/d^ 
are real and distinct, since generalization to the case of conjugate complex 

roots turns out to be direct^ (3) a = 1 in order to avoid its nuisance 

2 
value in the discussion • With these assumptions it is possible to have 

p « n in (3-23) with the denominator of each W (s) being a single monic 

factor in z^ that is, 

~sT 
1 ♦ c, e 

where c- may or may not be zero, but d^ / 0« There are n factors of the 
type (3-2ii) each representing a simple regression equation. In m of 
the factors c. / 0, in the other n-m factors c. « 0. The data storage 
associated with each ¥, (s) equals 2 when c / 0, and 1 when c. ■ 0."^ 
However, the input data that must be stored when c, / 0, is also the output 
data that had to be stored for the preceding cascade program ¥ ,(s)^ 
consequently, there is only one data to be stored for each W, (s) regardless 
of the value of c, , except for the first one W (s)o But when m <n (proper 
rational fraction), one e-^ - 0, say c^ = 0, making the total required data 



1 

If m^n, W(s) can be written as the sum of a polynomial and a proper 

rational fraction in z = e**^"**** The program corresponding to the polynomial 

part is a simple linear combination of input values* Discussion of this 

case is omitted without any serious loss of generality* 

If a^ / 1, only a simple multiplication has to be added to the program* 

Each W, (s) is the transfer function of just a regression equation and 

its data storage is the sum of degree of numerator and denominator, as 
discussed in the previous article « 
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storage n} a material reduction over the m -«' n data needed in direct regression 
programming. 

In order to translate the cascade scheme into an actual program, 
we may proceed as follows* First, we write (3-23) (with p = n) in terms 
of input and output transforms, as 



-OCs) .•o.Cs) "olCs) O'(s) 
T(s) 1^(s) 'i^(s) TCs) 

One way of making (3-2^) an identity is by letting 

l^(s) - 'l(s) 
I^is) « 0^(s) 

I^(s) = 'O^(f) 



^nf^) - Vlf^> 



which ma,ke 


0(s) = 


■ V^^' 






Using the relations 


(3-2^) and 


(3-26) 


we obtain 


OCs) 


"iCs) 


s;(s) 

'o^(s) 


o ft « 


TCs) 


its) 


^.if=> 



(3-2^) 



(3-26) 



(3-27) 



The various factors equal the respective program transfer functions^ namely. 
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-ii6. 




\(s) 


W^(s) - 


1 


A 


tCs) 


1 ♦ d^e-^^ 




cTjCs) 


W^Cs) - 


1 ♦ o^e-^^ 


\ 


^Cs) 


1 * d^e-^"" 




* 
• 






0(s) 


• 

' W (s) = 

n 


1 ♦ c e-=^ 
n 

1 * d e ^^ 
n 


J 


Vi(^) 



(3-28) 



where m of the c, are not zeno. Multiplying by the denominators changes 
the set (3-28) into 



-sT> 



(1 ♦ dj^e"^") 0^(s) =. I(s) 

(1 + d^e-^"^) -S^Cs) - (1 + C2e-"'^) o'^(s) 



-sTx - 



-sT. 



(1 + d3e-^') O3CS) - (1 ♦ c^e-^') ^^Cs) 



-sT, 



•sT, 



(3-29) 



(1 ^ d^e"^') -Ots) - (1 + c^e-^') ^.^(s) / 

The inverse transform of the foregoing set, -with one term of each equation 
transposed to the right side, is the desired set of regression equations. 

^ d^o-^Ct-T) 



•o^(t) « i(t) 



02(t) = r^Ct) 4' OgO^tt-T) - d^S^Ct-T) 
ol(t) * o^(t) ♦ C3O5 (t-T) - d^olCt-T) 



^t) 



« o; ^ct) + c^s: ^ci»T) - d;o(t-T) 

n-'X n n-l n 



(3-30) 
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The detailed coded program corresponding to (3-30) is shown in Appendix C» 

Cascade programming, although not referred to by that name, is 
a familiar technique in numerical procedures. However, the clear-cut and 
general equivalence of the direct regression and cascade programming is not 
always well understood • Cascade programming arises naturally from the kind 
of thinking prevelant in numerical work. Consider the simple example of 
solving the second-order differential equation, 

^ + P7 « a (3-31) 

dt"" 

The derivatives may be considered as the separate variables, y* (t) and 

y"(t)| then we obtain the following three sampled ftinctionss 

y"(t) - -?y(t - T) 

T'Ct) - Ty"(t) + y'(t - T) \ (3-32) 

'yCt) « 1Y»(t) +7(t - T) 

where the first equation of the set is derived from (3-31) while the second 

and third are elementary first-difference extrapolations. The set (3-32) 

indicates cascade programming because the output of the first equation is 

in the input of the second, and the output of the second equation is the 

input to the third. The peculiar thing in this ease is that the input 

to the first equation is not an independent function but directly related 

to the output of the last equation. This feature establishes the constraint 

imposed by the differential equation. 

The Laplace transform of the set (3-33) is 



IH = - e"sT^ 



°sT~ 



H = TY" ♦ e'^'-Y [ (3-33) 
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from which the e3q)licit relations between inputs and outputs are obtained 
as follows: 

Yti . - e "^^Y 

Y» « 11 Y" 

1 - e 
For realization by three cascaded factors, we have 

¥^(s) = -pe""^^ 



W,(s) - 



X 



1 - e"^' (3-35) 



-L "** 6 

It is clear that a single transfer function can be made to replace the 
cascaded system of threes thus 

¥(s) « W^(s)¥2(s)W^Cs) 



¥(s) 



-pT^e"*^^ (3-36) 



1 - ae"^*^ -t e~^®^ 



The corresponding regression equation is simply obtained as 

T(s) = (2 - pT2)e"^'¥(s) = e-^^^(8). (3-37) 

The inverse transform of this equation is 

Kt) * (2 - pT^)y(t - T) -^(t - 2T). (3-38) 

which could have been otained from (3-32) by the elimination of y*(t) and 

y"(t), but even in this simple case the process of elimination in the time 

domain is not direct. 
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The fact is that in niamerical work a cascade method such as (3-ii2) is 
much more generally used than the direct regression of (3-38). Often there 
is good justification for this preference^ for instance the values of the 
first and second derivatives may also be needed. However, when such or 
similar justifications do not exist, the direct regression may turn out to 
be simpler than cascading. In the present example, (3-32) calls for one 
more constant, two more multiplications and one more addition than (3-38). 
If the first two equations of the set (3-32) are combined, one multiplication 
is saved; furthermore, the manipulations in the direct method happen to be 
more awkward o Because in this case the input and output are the same quantity 
the formulas of Appendices B and C are not directly applicable, the requirements 
of the twD methods must be determined by actual trials. 
3«3i4. Parallel Programming 

If the transfer function of a program is e3q)anded by partial 
fractions in terms of z, (3-l8) tates the form 

f. fr, f 

WCs}-— i— -r- + ^ ^ni * .oo * — 2_^ (3-39) 

1 + d^e"^^ 1 4 d^e"^^^ 1 ♦ d^e"^^ 

as long as m <^ n. Thus, the transfer function ¥(s) is replaced by the 
stim of a nunfoer of simpler transfer functions| namely, 

W(s) « ¥^(s) ^ W^Cs) + ... + ¥p(s), {3*iiO) 

where some of the W, Cs) may be the combination of several partial fractions, 
but all are of lower degree than ¥(&) itself* 

The form of C3-40) may remind one of parallel combinations of 
network admittances. Paralleling means that the same input (driving voltage) 
is applied to all component admittances and the output (driving-point current) 
is obtained as the sum of individual outputs (current through each admittance). 
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The same interpretation can be applied to parallel programming • The programming 
will involve p regression equations all using the same input values, and all 
their outputs adding to produce the desired over-all output • 

To arrive at a more specific interpretation, we first make a few 
restrictions agains (l) W(s) is a proper fraction, i«e<>, m <^ n, and (2) 
the roots of the denominator polynomial are real and distinct* Then all 
constants f, and d, of (3-39) are real and in (3-UO) p can equal n| moreover, 
each term of (3-39) is a simple regression equation involving two constants 
and one data storage « Thus, the total number of data to be stored is only n* 
Just like in the case of cascade programming, the lower requirement for data 
storage of parallel programming may be a great advantage over the direct 
programming method* However, this feature does not mean that parallel or 
cascade programming should always be en^loyed in preference to direct 
programming. For instance, there is the case when m = 0| i«e*, the numerator 
of W(s) is 1 (or a )* Of the input values the program uses only the present 
one and the total data storage is n regardless of the programming scheme 
used| on the other hand, the number of constants will be n for the direct 
and cascade method, but 2n for the parallel method, putting the latter at 
a disadvantage. Similarly, if the denominator of the over-all transfer 
function lacks several terras (say, the denominator is 1 - b e ), then 
factorization of the denominator introduces all tenas, making the cascade 
and parallel program much longer than the direct program. Another factor 
which may militate against the use of parallel programming is the presence 
of multiple roots in the denominator* If a root is of multiplicity r, * 
it may produce up to r terms of degrees r, (r=»l), o»o2,l (the r-degree 
term never being absent) in the partial fraction exapnsion, but the same 
root will require only one r«degree, or r first-degree, cascaded factors. 
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In order to interpret the parallel method of prograTnining, we 
proceed in the usual manner. For the various terms of (3-itO) with p » n, 
we write 



W^{s) 


ICs) 


^1 


A 


1 ♦ djo-^^ 




WgCs) 


X(s) 

• 


^2 




1 * d^e-^T 


V 


W„(s) 


• 

0„(s) 


f 
n 




I(s) 


1 ♦ d e-^'f 
n 


J 



(3-ia) 



and 



¥(s) 



I(s) 

r(s) 



(3-1*2) 



Cross-multiplication by the denominators in (3-itl) yields the set 
(1 * d^e-^^)O^(s) = f^l(s) 
(1 + d^e'^'^)0^(s) "f^lCs) 



(3-k3) 



»sT, 



Cl * d e"""^)0 (s) « f rcs) J 

n n n — • 

while in view of (3-iil) and (3-li2)^ (3-iiO) can be written as 



0{s) « 0-Cs) + 0l(s) + ,,„ +0 (s) 
1 c n 



(3-Wt) 
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The inverse transforms of (3~U3) and (3-Uit) yield the desired set of 
regression equations, which follows o 

B-^Ct) - f^rCt) - d^3^(t-T) 
-OgCt) . fgiCt) - d^'S'gCt-T) 



0- (t) - f rCt) . d^S- (t-T) 
n n n n 

^(t) -^^Ct) ♦ -ouCt) ♦ ••. ♦ 'o (t) 
12 n 



(3-U5) 



The detailed coded program corresponding to (3-U5) is shown in Appendix D. 

Parallel programming has not been generally used in niimerical 
work. To the knowledge of the writer, the usual methods of niimerical 
analysis do not naturally lead from a direct regression equation, which 
has reference to several previous input and output values, to a set of 
simpler regression equations, each of which refers only to the last 
input value and to a preceding output value* By the method of frequency 
transformation the parallel method is found quite directly. 



In case of complex d, *s in (3-39) ^ a combination of two conjugate complex 

partial fractions in z will resxilt in a slightly more complicated regression 
equation, involving one additional input and output. 
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3 '35 Compaii-son of PrograMiiing Methods 

The purpose of this section is to compare the effectiveness of 
the various methods of program realizations based on the transfer 
functions of the programs* A complete general treatment appears too 
far-fetched and, therefore, this study is limited to a certain class of 
programs. Despite these limitations, which are discussed below, the 
investigation is sufficiently general to show how the results can be 
used to improve the instruction code of a general -purpose computer or 
to design a special -purpose computer, when these are used in control 
applications* 

The three methods which will be compared are listed below: 

(a) direct programming, 

(b) cascade programming, 

(c) parallel programming. 

Other programming schemes may be derived from the rational transfer 
function W(s). One may carry out the long division in z of the numerator 
by the denominator until he arrives at a certain number of terms of the 
quotient. The transfer function can then be expressed as the sum of the 
quotient terms and of the remainder divided by the divisor (the original 
denominator). Any number of variations can be obtained by stopping the 
long division after different number of steps, but only in the most iinusual 
cases can this approach be expected to yield a more efficient scheme of 
programming than the three major methods discussed in the preceding sections* 
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Other schemes that are even more artificial than the long-division scheme may 
be derived, but no other general programming method has been fo\ind that ^ves 
promise of effectiveness conqjarable to the three which are considered* It 
is noted that in certain cases a combination of two of the three listed 
methods may turn out to be more efficient than any one« An example of 
such a case is described belowo 

As the basis of comparison of programming methods, the requirements 
in storage and time are used* The particular application or purpose ctecides 
which of these two factors should deserve morei attention* It is assumed 
that the complete sequence of instructions, as used at each sampling point, 
is store d| the possibility of cycling programs, which re-uses a short 
sequence of instructions for the calculation of each term, is not discussed* 
Essentially the Whirlwind I code is used throughout, but variations are 
considered* 

As a starting point we recall that the transfer function of a 

linear program is, 

•»sT ^ -2ST ^ ^ -msT 
Tsf \ a ^^ a_ e + a^e ♦•••♦e.e /oi/iN 

I(s) 1 ♦ hjB # b^e + .•* ♦ b e 

In general, m and n may be any positive integeiss and indeed, their relative 
sizes will hardly influence the comparisons to follow* Nevertheless, it is 
helpful to distinguish three cases? 

(l) n » Go (3-1^6) reduces to a polynomial in z - e j i*e*, 
the new output value depeads oaly on present and 
past input values, not on past outputs also* 
Present-day numerical analysis abounds in numerical 



1 

This is in contrast with networks where certain restrictions on the 

degrees of numerator and denominator pol3momials often exist* 
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processes correspondir^g to this special case* 
(2) m < n«. (3-1*6) has the form of a proper rational function 

of z in this case. In Sections 3 032, 3*33^ and 3«3i* 
dealing with the various prograimiing schemes, this 
case was assumed for the sake of simplicity* 
C3) m^n. The rational function in z of {3-k6) may be called 
improper, but it can be converted to the sum of a 
polynomial (Base l) and of a proper rational fraction 
(Case 2) in z =» e * 
In order that the storage and time estimates to be arrived at should apply 
to all cases, it is necessary to define the folldwiijg quantities m.th 
reference to (3-ii6)j 

ffi « degree of numerator, 
n « degree of denominator. 






= one less than, the number of non-zero constants in 
the n\imerator (m^ ^ n) • 

= one less than the number of non-zero constants in 
the denominator (n. ^ n) • 



m = one more than the excess of m over n| i.e#, 

m^ = m-n*l when m ^ n, and m « otherwise* For 

proper rational fractions m < n and m « OC 

e 

On basis of the coded pr*ograms shown in the appendices, the table of 
Fig, 3»U summarizes the storage and time requirements in terms of the 
quantities just defined* This tabulation is more general than the results 
given in the appendices, for in the appendices it was also assumed that 
none of the constants were zero, that is, m^ = m and n^ = n| furthermore. 



1 

Examples are numerical methods based on polynomial approximations with 

equidistant spacing of the independent variable « Indeed, such examples 

form not an insignificant portion of the available numerical techniques* 
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only case (2) was treated making m = G* On the other hand, in the tabiilation 
of Fig« ^•h these restrictions of the appendices are absent, but the following 
assxmptions are still made: the roots of the numerator and denominator 
are real and distinct, and the straightforward prograraming techniques of 
the appendices is used» Thus, the constants stored are those that appear 
explicitly in the various regression equations. Actually some saving in 
instructions would result from the use of certain ratios of these constants. 
For instance, the regression equation fcf* (^^^h^Tj 

'S^Ct) « f-jT(t) - d^o^(t-T) (3-1*7) 

takes six instructions, as shown in the coded program of Appendix D, 
If, however, C3-i*7) is written as 



S^Ct) - f^ 



t(t) - -5^ %(t-T)J , (3-U8) 



its coding would cost five instructions only, but certain questions on the 
relative sizes of the constants would arise. It seemed best to avoid such 
questions, because the considerations here are rather general and the value 
of a too-specialized treatment is questionable. 

The comparison of the three methods of prograraming can be undertaken 
by considering each item of Fig, 3»lt» Because of the straight sequential 
programming the time requirements are the same as the storage for instructions 
and, therefore, consideration of storage will give a complete picture. 

As far as the number of constants stored are concerned, the direct 
method is not worse than the cascade, which in turn is not worse than the 
parallel method. This is so because in the direct method only the non-zero 
constants of (3-ii6) have to be stored, while factorization in the cascade 
case will produce as many constants as there are roots in z. In the parallel 
method two constants (root and residue) are produced for each denominator 
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root and if the numerator is not of lesser degree than the denominator, 
further terms and constants result. As an example, consider 



¥(s) 



^ 1 ^-sT 

(3-U9) 



1 . 3 ,-2sT ^ 1 ^-UsT 



for which 



m - 1, \ " ^» Q 



m « 



n » U, n, ■ 2, 

According to the table of Fig« 3»U the various constant storage requirements 
are 

directs m. ♦ n. + 1 = I4 

cascades m + n ♦1 = 6 

parallels 2n = 8 
These figures can be simply checked. In the direct case the four constants 
are apparent in f3*^)«» For the cascade case, the transfer function is 

written as 

$ 1 ^-sT 
T 1 1 Ti " T 

"^=^ = , , 1 -sT • :; — T-zsr - ., , 1 , -st • :; — r-=if — » f3-5o) 

l*;z.6 i-_e 1 + — e 1 e 

n \l2 \2 2 

and the six constants in question ares ♦l//^, -l/lT^, *l/2, -l/2, ♦^/S, 
and -l/iio The manner of programming illustrated in Appendix G actually 
necessitates the separate storing of positive and negative constants, 
even though of the sam6 magnitude • 

For parallel programming ¥(s) of (3-it9) is expanded in partial 
fractions in terms of z and takes the form 
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$ ♦ 2^2 


^ - 2i/r- 


W(s) = 


8 


, 8 


^ 'k-' 


1 1 «-sT 




9 


1 




l^fe-^ 


1 2 ' 



(3-51) 



The eight constants to be stored are evident in the foregoing equation. 

The next item of comparison is the data storage ^ for idiich the 
above example reads, on basis of Fig. 3<,J4 

directs m ♦ n = 5 

cascades n « I4 

parallels n ■ U 
These figures can be verified in the three foregoing equations. The numerator 
of C3-I49) indicates that one past input value (corresponding to the e*^ term) 
must be stored^ the present input is used as it arrives and then stored 
as the past input for the next calculation, as shown in Appendix B. Thus 

the numerator implies one data register only. Similarly the denominator 

-sT 

implies the storage of four past output values, even though the e and 

«=3sT 
e terms are absent^ for the corresponding past outputs must be remembered 

for the next calculation. 

For the cascade method, (3='5»0) seems to indicate five past data 

>-ST y- ■ 

to be remenberedi however, the e term of the last numerator rehire id 
a past input that is also the past output of the preceding factor, since 
in cascade programming the input of a component program is the output of 
the previous one* 

In case of parallel programs the four past data are quickly 
identified with the e terms of (3-5l)ff 
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The expressions for instruction storage and for time requirements 
are identical, and produce the following tally in the present example? 

directs 2(ra ♦m,'»'n-»'n,) + 7»23 

cascades 3ni ♦ Un '^ 6 ■ 2^ 

parallels 7n * U » 32 

No verification of these figures is carried out by detailed coding of the 
programs because the appendices cover the general case* The advantage 
seems to be on the side of direct programming as far as time is concerned, 
but this advantage is slight and arises from the fact that in the present 
ex&mple two denominator constants are zero. An advantage of direct programming 
appears also in the total storage requirements for the same reason: 

directs ii ♦ ^ ♦ 23 - 32 

cascades 6 ♦ H ♦ 25 ♦ 1 = 36 

parallels 8-»'it*32 4'l«i45 
This example, as well as the tabulation of Fig. 3 ^I^^ 'indicates 
the disadvantage of parallel programming. It seems that this kind of 
programming may have an advantage over either of the other two in certain 
cases, but hardly ever over both at the same timeo Thus, the choice narrows 
down to direct and cascade programs, or possible combinations thereof. 
To show how a combination of methods may be used, we write (3-ii9) as 

5,1 g-sT 

W(s) « ,^ .m . ^ " ^ ^ .,« (3-52) 

1 «= I e"2^^ 1 » I e'^^^ 

which indicates a cascade combination of two direct programs, for which 

respectively 
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The direct program of each cascaded coiRponent is somewhat simpler than it 
would be for two separate direct programs because the input and output devices 
are manipulated only once for the composite program, rather than once for 
each component programo This saving amounts to six instructions, thus the 
instruction storage or time requirement iss 

first compqnentss 2(m + m. + n -^ n, ) ,+ 6 = 12 

second components! 2(m ♦ m. + n + il) •♦• 6 « 16 

saving as indicated above -6 

total instructions '22 



Four constants appear in (3-52), two of which are accidentally identical, 
and one of which is made 1| thus, the constant storage iss 
first components m, ♦ n^ ♦ 1 « 2 

saving » -1 1 
second components m, ♦ m, ♦ 1 ■ 3 

SAving « ^ 2^ 

total constants 3 

A saving arises in data memory also, because the past input of the second 
component is also the past output of the first one» This gives the following 
need of data storages 
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first component § 


m ♦ n » 3 


second components 


m + n = 3 


saving 


-1 




h 







The results of this example are summarized in Fig* 3o5^ which shows a small 
advantage of the mixed method over the direct one« 

To pursue further the detailed comparison of these various methods 
of programming would lead to undue specializations in the Whirlwind code and 
to results of doubtful general value o The illustrated attack on the 
realization problem^ however, shows how a useful estimate of the complexity 
of coded programs can be gained from the evident properties of their transfer 
functions* Three further problems will be touched on briefly? (l) computing 
delays^ (2) means of using the results to select computer codes^ and (3) 
means of using the results to design special ^purpose computers. 

A consideration that has been omitted in our discussion is the 
delay incurred through the computation itself o If a digital computer is 
used as part of a number of control systems =- say, ^0 systems — , then 
in each sampling interval it performs ^0 computations, one for each system. 
The time of a computation is then at most l/^O of the sampling time, T, 
and this delay is presumably negligible « If, however, one digital computer 
were used with each system, the computation may and, for the sake of efficiency, 
should take an appreciable part of the sampling time. Such a delay would 
be very serious and the computer would have to perform a prediction in 
addition to the required compensation. In turn, this would lengthen 
the program, make it less effective^ and may even force a longer sampling 
time| indeed, in a marginal case, in which the original compensating program 
had a delay nearly as large as the sampling time, the effect may become 
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cumnlative, since a longer sampling time would in turn require a better 
and longer program, and so ono In such mai^inal cases and in any case in 
which the computing time is not negligible with respect to the sampling 
time, the direct programming has a tremendous advantage over all other 
methods o A glance at the direct regression equation (3-20) shows clearly 
that ail terms but the first one on the right side of the equation can 
be computed before the new input value is obtained* The computing delay 
will thus be the time of merely calculating the term, a i(t), and adding 
it to the already prepared partial result <» This delay may conceivably 
be negligible • 

All realizations of real-time linear programs involve accumulation 
of products as their arithmetic action and the transfer of data from one 
register to another as their manipulative action. In case of a single-address' 
instruction code, such as that of ^feirlwind, the ex (exchange) operation 
was shown in Appendix B to be very helpful in improving the efficiency 
of the code© Other in5>rovements aj^e possible by incorporating special 
©peratibns which facilitate the particular type of programs on hand* 
Computers using multiple =address codes could be particularly efficient in 
such applications o For instance, in a three-address code an instruction 
could locate a constant, a piece of data, and transfer that data to a 
third address, after which it would multiply the constant and data 



1 

The second composite program in Appendix B is written in this manner* 

2 

Each instruction specifies an operation and the storage address of a 

single operandi 

This operation exchanges the contents of the accumulator register with 
the specified storage register* Thus, one instruction performs a double 
duty by obtaining aew data from storage and ^so transferring to storage 
a partial result* 
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acctunulating this product with the partial resiilt always left in the arithmetic 
element of the computer o This single order would complete both the arithmetic 
action (accumulation of products) and the manipulative action (transfer of 
data to an "older-data" register) associated with one terra of a regression 
equation o 

Similar considerations allow one to adapt special -purpose or 
fixed-program digital computers to control specifications. To be somewhat 
specific we assume that the conputer is used as part of a single control 
system and will have to perform only one computation in each san55ling period o 
The coi!?)uter would not operate appreciably faster than one computation per 
sampling period and in order to minimize the computational delay it would 
follow a direct regression program. In order to keep such a single-system 
computing equipment from becoming excessive, a serial computer woTild 
probably be usedo The program of the computer would be fixed to correspond 
to a direct regression program of certain complexity as defined by the 
degrees of m of the numerator and n of the denominator of the program 

transfer function « The constants could be set manually on toggle switches 

2 
or relays^ or they could be stored on the same high-speed storage device 

on which the data are stored. A serial adding unit with proper switching 

equipment would allow the multiplication of constant and data (by repeated 

additions) and the addition of such product to the accumulated partial 

result. The physical size of such a digital control tinit may be quite 

feasible in certain applications and the design of such a siuple special-puipose 

digital computer would be particularly ^^stified if the incoming data were 

sampled and digital to start with. 



1 

A serial computer operates on each digit of a number in sequence^ thtis^ 

the equipment is not duplicated for each digit. 

2 
Magnetic-drum memory, for instance. 
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3«'^ Synthesis of Programs in the Frequency Domain 
3oitl General Ssmthesis Procedure 

The synthesis of computer programs in the frequency domain may 
be broken down into the three following stages Cl) specification of the 
desired frequency characteristic or lous of WCjoi), (2) approximation of 
the desired locus by a realizable program transfer function, and (3) 
realization of the program. One way to determine the desired locus is 
from the Laplace transform of the operation the computer is to perform* 

The second step is the difficult part of the problem. The desired frequency 

-sT 

characteristic must be approximated by a rational function of e .No 

general rules are available for making this approximation, but before 
making the approximation, one should gain some experience in analyiZing 
program building blocks in the con^jlex plane* Possibly the most systematic 
approach, at present, to the approximabion problem is to make successive 
approximations to the desired characteristic, using the basic program 
building blocks of Section 3»2« The third step involves only a straight- 
forward inverse Laplace transform. As an example of program syhthesis 
in the frequency domain, a program for differentiation will now be synthesized. 
3oi^^ Synthesis of a Differentiation Program 

An ideal differentiator establishes the following relation between 
input and outputs 

OCt) = 1^ ^^^K (3-?2) 

Disregarding initial conditions, the Laplace transform of (3-52) is 

H(s) . ofci = s, D-53) 

Ks) 

and this is the desired transfer function* For s « joo, H(s) becomes 

H(ja) - j6). C3-5il) 



Report R-225 



-65- 



So the locus of the desired transfer function is the iraaginaiy axis. This 

completes the first step of the synthesis procedure* 

-sT 
vThe second step is to find a rational function of e that 

approxiinates this locus. This approximation is to be made by geometric 

considerations based on the desired locus o In this particular exanple it 

is also possible to employ analytic considerations based on the desired 

transfer function of (3-5U)« It so happens that in the present case the 

analytic approach is simpler^ nevertheless, the geometric approach is shown 

first. 

The crudest numerical approximation to a first derivative 

is the first divided difference. 



^(t) , iCt) - TCt-T) 
T 



-sT 



The Laplace transform of (3-55) is 

T 
Thus the transfer function of the differencing process is 

^t e ^ Sts) ^ 1 - e"®^ 
^ TCs) 



(3-55) 



(3-56) 



T 



(3-57) 



Figo 3*6 shows the locus of W (jco) and compares it with the desired one. 

At low values of ©T (i;ei, -idien the frequency of the input ftmction is Xow 

''^i(t) -i(t-T) 



s "locus 




locus 




> * 



Figure ^•6 Comparison of first derivative and first difference operators 
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with respect to the sairipling frequency) the two loci agree reasonably well* 
If we could straighten out the circular locus, we would have a better 
approximation of the desired locus. Figure 3»7 illustrates a geometric 
constniction that straightens out the locus and gives us ideal phase 
characteristics » 



-sT 



di(t-»T) 
dt 



i(t) - i(t-T) 




(a) Frequency domain 



(b) Time domain 



Figure 3*7 Derivation of an ideal phase> realizable differentiation 
operator 
The vectors (l/T){l - e*""^®^) and (l/2)Cl + e"'^® ) are drawn for 
a particular frequency* Using the geometric rule that a triangle inscribed 
in a semicircle is a right triangle, one can readily show that /o^/ * \^\ 
add up to 90 « However, since c< is a positive phase angle, it must be 
subtracted from p (which is negative) to give a resulting angle of -90 ^ 
isghich is the phase of an ideal integrator* It follows that division of 
the P^locus by the o<. -locus will yield an ideal-phase formula* The 
resulting transfer function is 



W^Cs) « I i 



-sT 



1 ♦ e 



^ 



{3-^8) 
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which has the desired phase in the range, ^ooT ^ir» The interpretation 

in the time domain is both plausible and illuminating. The inverse transform 

of (3*='58) shows that, 

-yCt) * y(t - T) . i(t) - Kt-T) ^^^^^^ 

(3'°^9) states that the average of the derivatives at two neighboring 
points is approximately equal to the divided difference for those points. 
It is interesting to note that the same approximate transfer 

function, (3=58), can be obtained analytically based on a rather good 

„sT 

approximation for e • 

^.sT ^ 2 (3^Q) 

l.isT 

Solving (3=60) for s yields 

3^4 1 -" ® , , (3-61) 

•L ^ e 

Although in this particular case the above analytic approach is 

simple and fairly accurate, its general use has certain drawbacks. The 

most obvious one is that the rational function of "s" to be approximated, 

which in the present case is "s" itself, is in many cases not explicitly 

knowni rather it may be obtained as an ^proximation to a desired locus 

or ai?55litude and phase response. Then to approximate the rational 

function of *'s", ^ich itself is but an approximation, by a rational 

"sT 
function of e puts the designer on shaky grounds, and it might lead to 

far more involved programs than necessary. There is no substitute to 

going back to the original specifications and designing directly on their 

basis. Another disadvantage of the above analytic approach is that it 

is not general. One could replace all **s" by the approximation (3-6l), 

but how one would get a better solution is not obvious. 
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We have an approximation of the differentiation operator, so the 
next thing to do is see how good it iso Since the desired locus and its 
approximation lie along the same path, a locus study does not give a good 
comparison. In such a case separate amplitude and phase plots can be studied* 
For s « jo), WpCs) becomes 



WgCjco) = ^ 



.coT 
2 1 - e'^"'* _ 2 e'*'^ 2 



1 ♦ e 



gosT 



■ e^T 



.coT T ^ 2 ^ 



(3-62) 



which verifies the previous statement that WpCjo)) has ideal phase characteristics. 
Hence, it is suff icent to study the amplitude characteristic only. 

H(jcd) - 0®^ (3-63) 

therefore. 



W2(J®) 
H(jco) 



tan- 



o)T 



2 



(3-6U) 



Thus we see from (3-61i) that the ratio of the approximate function to the ideal 
one is always greater than unity. Figure 3 08, a plot of the amplitude 
characteristics, shows us that for low values of gjT, say for coT ^—tk 9 the 

A 

I (1) H(j6)) = 03T 

coT 



(2) W^Cjco) « 2 tan 1^ 



(3) W^Cjo)) « l.eiii tan ^ 




asT 



ir o)T 

Figure 3 -S Comparison of amplitude characteristics of Differentiating Operators 
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differentiation program will give quite good accruacy. For certain control 

applications, values of <aT up to -z or even ^ might give acceptable accuracies • 

An examination of the amplitude characteristic of W^(s) in Fig» 3»8 

reveals that if W^Cs) is multiplied by a constant, which is slightly less than 

unit, we will obtain a better de±*ivative on the average. The new transfer 

function is 

-sT 
W^Cs) - C W^Cs) « C I ^^^=-~~ip • (3-6^) 

X T e 

tr 



Let us arbitrarily choose C so that ¥^(s) - H(s) for sT = j -r* Then, 

tr If , 



2 C tan I « I I (3-66) 



so 

tr 



6 tan "T 



*2p- - 0.907* (3-67) 



The improved transfer function is 



W (s) « 1^ ^ ' ^'% s (3-68) 

and its amplitude characteristic is also shown in Fig, 3 ,80 

Both curves 2 and 3 of Fig. 3.8 accentuate high frequencies which 

may be present at the input because of noise <t In this case, a transfer function 

whose amplitude characteristic is like that of curve k would be a more desirable 

approximation for differentiation. 

The inverse trsoisform of (3-68) completes the synthesis of a 

differentiation program. The result is 

S(t) - 3^ (Tct) - i(t-T)J - cTCt-T). (3-69) 

The accuracy of this differentiation program may be determined from Fig. 3*8, 

Curve 3« 
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CHAPTER IV 
FREQUENCY ANALYSIS OF SOME NUMERICAL INTEGRATION FORMULAS 

In this chapter we apply the methods of frequency analysis to 
several nxuaerical integration formulas: the trapezoidal, Sin5)son's 1/3 
rule, Simpson »s 3/8 rule, and Weddle«s rule* Frequency analysis is 
applied to determine the stability of these formulas, compare their ac- 
curacy, and compare their transfer functions with that of the ideal 
integrator. 

kol Numerical Integration 

In the numerical integration of definite integrals, the range 
of integration is divided into a convenient number of equal intervals, 
and the values of the integrand are defined only at the ends of these inter* 
valso Essentially this is the same as sampling (or impulse modulating) 
the integrand* Let the distance between sair5)les be T« To obtain an 
approxiuiate value of the integral we may determine an n order polynomial 
that passes through n + 1 of the sampled points and integrate the poly- 
nomial over the corresponding range, repeating the process until the com- 
plete range of the original integral has been covered* If the sampled 
points are joined by straight lines, (approximation lay a first order of 
poljmomial) the resulting integration formula is known as the trapezoidal 
rule (each interval of the integrand is approximated by a trapezoid)* 
Joining the points in each group of three sampled points by a parabola 
leads to an integration formula known as SiBf5)son»s 1/3 rule* If the 
points in each group of four sampled points are joined by a cubic curve, we 
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get Simpson* s 3/8 rule, •phe trapezoidal rule and Sin^son's 1/3 rule 

are quite widely known and used, but there is another one, called Weddle*s 

rule, that is used to obtain great accuracy. Joining the points in each group 

of seven sampled points by a sixth order polynomial leads to Weddle's rule. 

In each case the range of integration should be divided into an integral 

multiple of "n** intervals. For example, to use Weddle's rule, the range 

of integration should be divided into 6, 12, 18.. ,,,. equal intervals. 

In i^iat follows we shall designate the transfer function of an 
ideal integrator as H(s)s Thus, 



H(s) = i 



(k-1) 



with which the approximate integration formulas will be compared. 
I|..ll Trapezoidal Riale 



Using the trapezoidal rule the definite integral. 



o(t) = 

may be approximated 1:^, 



1 



i(x) dx <? 



(U-2) 



oo 




o^(t) ^^/ / i(t) + i(t - T) 



]. 



i(t » T) + i(t " 2T) 




The Laplace transfona of (h-3) is 
0,Cs) = | [ 



(1+ e"^^) (1+ e'^*^) + e'^^^-j. ) 1 l(s). 



(h-k) 
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So the transfer function is 



W^Cs) = -fey = I i^ (W) 

A little algebra shows that for s = ^^ 

-HOST " ^ """^ r* ^^-^^ 

l4»12 Simpson ys 1/3 Role 

Using SiBipson's 1/3 rule the definite integral (U-2) may be 
approxiiaated by 

OgCt) = I J [i(t) + lii(t - T) + i(t - 2T) J + 

[i(t - 2T) + i|i(t - 3T) + i(t - kT) J+. ...... i (U^) 

•flie Laplace transform of {I4-8) is 

OgCs) =1 I(s) (1+ h e-'^ * e -^^^) (1+ e-2^f + e -^^%... j 

(U-8) 



or< 



"•sT "»2sT 

1 - e 
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lEherefore^ the transfer function for Simpson* s 1/3 rule is 
OpCs) T . -sT -2sT 

^2(-> = 4riT = J ' ,!^-/st' — • (^-^°) 



Dividing (Ji-10) by (I4-I), letting s = j<a and using some algebraic and 
trigonozaetric loanipulations leads to the ratio 

^2^^^^ _ caT 2 -¥ cos eaT n ^^^ 

"hO^ - T sin<»T • ^^-^^^ 

I1.0I3 Simpson's 3/8 Rule 

The approximation toihe definite integral (l(-*2) that is obtained 
msing Simpson's 3/8 rule is 

03Ct) = ^J [ i(t) + 3i{t - T) + 3i(t - 2T) + i(t - 3T) J + 

|i(t - 3T) + 3i(t - IjT) + 3iCt - 5T) + i(t - 6t) + J 

(1^-12) 
The Laplace transform of ©oCt) is 



O3CS) = |2 i(s) 



(1+ 3 e -«^ ^ 3 e-^^^ ^ e^^^^f^^ e'^^'' ^ e^'\,J 



(k-13) 
or^ 

A /'o\ - 3T 1 4. 3 e"^*^ + 3 e"^^^ -¥ e°^^^ ^f. .,t,. 

0^(s) = ^ ^^ ^^_ I(s) Cij»il^) 

X «• e 
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Hence f or Siiapson's 3/8 rule, the transfer function is 
V M^) -^T 14. -? e-s^+ -i e-2sT -SsT 



For s :s j«a, the ratio of W-{ jea) to H( joa) is 



W3(j(o) ^ 3 ^T 1 + cos qdT 

"^^^^ " "^ ' (1+ acoseoT) tan^ * 



(14-16) 



A considerable amo\mt of manipulation is reqiiired to obtain the above form. 

'« 

UclU Weddle^a-Rule 

Ely Weddle's rule the approximation of the definite integral 
(ii«2) is 

©l^(t) = ^ j [i(t) + ^i(t - T) + i(t - 2T) + 6i(t - 3T) + 

i(t - ItT) + JiCt - 5T) + i(t - 6t) J + k(t - 6T) + 
5i(t - 7T) + i{t - 8T) + 6i(t - 9T) + i(t - lOT) + 
^i(t - IIT) + i(t - 12T) 



] 3 



(lt-17) 
In the same manner as before, the transform of Oi (t) is 

Oj^(s) - j5 :55| i(s); 

1 - e 

(lt-18) 
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so the transfer function for Weddle»s rule is 

V^^ " "rTiT " 15 iTT^sf 

(U-19) 
Biy using a considerable amount of algebraic and trigonometric manipulation, 
we get for s = jcos 

Wj^(j(©) 3 ajT 1 4- 3 cos (aT ■» cos^ (oT ^ i «rt^ 

~~~~ (1+ 2 coscdT) sina>T • ^^"^°^ 



U»2 Coisparison of Numerical Integration Formulas 

With the above equations, we can get a complete picture of the 
four approximation formulas in both the time and frequency domains. 
Equations (it«5), (b-10), (l^-lj), and (U-19) are the transfer functions 
of each of the numerical integration processes and from these the stability 

of each one can be determined* Let us now examine the denominator of each 

-sT 
transfer function. If the change of variable, z = e , is made, it is 

easily seen tiaat the magnitude of the roots of all the denominator poly- 
nomials is unityi however, there are no multiple roots. Therefore, each 
of the numerical processes is stable. 

Now we must consider the accuracy of each of ihe integration 
formulas. Eq\iations (U-6), (1^-11), (1^-16), and (U-20) give the ratio of 
the particular transfer function to that of the ideal integrator. In Figure 
itol these ratios are plotted as functions of coT, and we see clearly that, 
of the four, Simpson* s l/3 rule and Weddle»s rule are the best for coT 
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Trapezoidal Rule 
Simpson* s l/3 Rule 
Simpson's 3/8 Rule 
Weddle's Rule 
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below ^ radians • For example, suppose that we wish to integrate a sine 
wave of radian frequency on and want the error to be less than 2.5^. 
For each of the approximation fommlas, how many samples laast be taken 
in a cycle of the sine wave? The answer can be obtained rapidly from 
Figure k»l by noting the frequencies at which the amplitude ratios become 
0«975 or 1»025, as listed below^ 

Trapezoidal oa T = 30 5 12 samples/cycle 

Simpson »s 1/3 Rule co^T = 7^ 5 5»8 » " 



o 



.0 



Simpson's 3/8 Rule 09 T = 60 j 6,0 



o 



Weddle«s Rule ffl T = 80" j k*S " 



tt 



o 



The number of samples per cycle is indicated for each rule, and this is 
obtained Tijy dividing 360^ by the indicated angle • 
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APPENDIX A 



Proof that the Locus of Q( ja>) Crosses the Real Axis either Normally 
or Tangentially at co » and^ 2 

*■ \ ~sir 

Recall that Q(s) is given by the polynomial in e , 



QCs) « > b, e'^'^ $ b « 1 . 

^ 7i k ' o 

k « 



For s = and - j ^ , Q(s) is real because each term of the polynomial is 
real. Since the locus is symmetrical about the real axis^ it must cross the 
real axis at these points • 

In order to examine the behavior of the locus of Q(jco) at these 
points, take the derivative of Q(s) with respect to s« 

-sT 



# '^ -k Tb^ e' 

ds X k 



dQ — sT 

Observei that -rr is also a polynomial in e $ therefore, it will also be 

real for s « 6 and — j "^ • 

Now consider the derivative in the neighborhood of s » 6 and 

♦ j^Q-. If j^ / at these points, we will prove that the Q-locus crosses 

i 
the real axis perpendicularly. In the region of interest let ds * j ^ , 

where cf is a small increment of (»• Since ^ must be real (and unequal 

to zero as we have assumed), dQ = j^j |dQ| in order to ^ktlaty this condition. 

(a<.E*D.) 



¥e must now discuss the case in \diich -sr " for s « or — j#" . 



d§! _ n ^ - -. n «« ±4J3- 

ds 



First; observe that if j- * 0, Q must have a saddle point in the region 
near the point where ^-^ « 0. Let us now make the change of variable. 



1 '■ 

For an excellent discussion of the l&ehavioui* of functions near saddle 

points, see Guillemin, "The Mathematics of Circuit Analysis^f JJbfiri'Mley 

and Sons, New York, 191^9, pp. 298-302* 
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'79= 



'sT 



^ so that Q(s) becomes D(z) which is 

111 



D(z) = J 






In the iramediate vicinity of a saddle point, the function behaves as 

D(z) ^ G •► C (z - z)^ 
^ o p 

in which the C's are constants, z is the value of z at which the saddle 

o 

point occurs, and "p - 1" is the order of the saddle point. In this case, 

z « ♦!« In plotting the locus of D(z), we map the \mit circle of the z-plane 



into the D-plane (see Figo A-l)o 
j 7 >^ z-plane 




Consider the map in the vicinity 

of a possible saddle point (z « ♦!). 

Observe that for z near a ; 

z - z « dz ♦j /dz I . There- 
fore, in the vicinity of a saddle 

point DCz) is 

D(z) = C^ ♦ C C+j |dz |)P. 



unit 
circle 



Figo l^l Unit circle in the z-plane that 
maps into the I)~plane 

This readily shows that if **p" is even, the locus in the D-plane (or Q-plane) 

is tangent to the real axis» If "p" is odd, the locus is normal to the real 

axis« 

We will now summarize the results obtained* 



dQ 



JX 



a) If 3^ f^ for s « or *j"y", the locus of Q(jci)) is normal to 



fl 



the real axis for s « G for ♦^•^ respectively 
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b) If ^ = for s « or + j-^ Q(s) has a saddle point at the 

-sT 
point where the derivative is zero* The change of variable, z « e , 

permits us to write, D(z) « C + C (z - z )^ for z in the immediate vicinity 

of the saddle point* If *♦?" is even, the C^plane locus is tangent to the 

real axis at the saddle point. If "p** is odd, the locus is normal to the 

real axis. 



Report R-225 

APPENDIX B 
Coding of Direct Regression Programs 

The regression formula 

^Ct) = a i(t) ♦ a^rCt-T) ♦ ... ♦ aT(t-mT) - b.o'(t-T) - ... - b 'oCt -nt) 
1 ml n 

(B-l) 

is to be prograinraed. Ass-ume that the data and the constants are stored 
as follows s 



Register 


Content 


Register 


Content 


No« 


(Constants) 


No. 


(Data) 


A.O 


a 




I.O^ 


Kt) 


A.l 


^ 


I.l 


T(t - T) 


• 


. 


• 


• 


• 


• 




• 


« 


• 


• 


• 


A.m 


a 
m 


X.m 


i:(t - mT) 


B.l 


-\ 


0.1 


Z(t - T) 


B.2 


-b2 


0.2 


^t - 2T) 


• 


• 
• 


« 




« 


• 


. 


. 


B.n 


-b 
n 


Om 


'oCt - nT) 






R.O^ 


Partial and 
final results 



These registers are not used in the second composite program. 
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The program will first be coded in two distinct parts 8 arithmetic and 
manipulative o The arithmetic part performs, at each sampling point, the 
arithmetic operations called for by the above regression formula and thus 
calculates a new output values 





First Program^ 


Arithmetic Portion 


Register 
No. 


Content 
(Instruction) 


Result 


P.l 
P*2 
P.3 


ca O.n 

mr B.n 
ts R.O 


4. -b ^(t - nT) 

^ n 


P«k 
P-^ 
Po6 
P.7 


ca 0*n-l 
mr B.,n-1 
ad R.G 

ts UnO 


-i^ -b ^"^ (t - (n-l)!fl -b 15Ct - nT) 






« 
« 
• 


• 
• 
• 


• 
• 


P-(im-ii) 
etc. 


ca 0.1 
mr B.l 
ad R.O 
ts E-0 


n 




k=l 


P.(l;n+3) 


ca X.m 
mr A.m 
ad R.O 
ts R.O 


n 
_V a i(t - mT) - N h.TSit - kT) 




f m ^v =■ "k 

k«l 



The code is explained in Sc-D, Thesis "Treatment of Digital Control Systems 
and Numerical Processes in the Frequency Domain," J.M, Salzer, Appendix 1.C, 
¥ol. 2, August 1, 1951, M.I.T. 
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Contimieds 



Register 
Noo 


Content 
(Instruction) 


Resxilt 


o 


« 






P.UCn-hn) 
etc© 


ca loO 
rar A«0 
ad ReO 
ts ^0 _ 


-^ o(t) 


Po(l4n+Um+5) 


si 

re R,0 


selects the relevant output device 
(as specified by the address section) 

records output, o(t), into output 
device 



It is clear that in this illustration each term of the regression equation 
costs h instructions* 

After the calculation of o*(t) at a particular saii5)ling point 
the data storage has to be rearranged for the next calculations the present 
"Sft - nT) can be lost, all other '£f(t - kT) are to be stepped dowi one storage 
register, and the new output value, ©"(t), just computed is put into 0«1| 
the rearrangement of the i(t - kT) is analogous, and the new input value to 
be^3?eceived goes into I«0o The coded program performing these manipulations 
folloi^Sc 
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FIRST PROGRAM^ MANIPULATIVE PORTION 



Register 
Noo 


Content 
(Instruction) 


Description 


P<,(Iun+lm+6) 


ca Oon-^l 
ts 0<.n 


/ moves 'S' t •=• Cn=l)T into location 
J of TCt - nT) and loses o(t - nT) 




ca Ocn<=2 
ts 0.n=l 


T moves 0^ ft = (n«=-2)T] into 
j o [t - (itl)'g location 


o 


« 
« 






ca Ool 
ts 0o2 


~) moves o(t « T) into 'o(t - 2T) 
J location 




ca RoO 
ts 0,1 


1 moves o(t) into o(t - T) 
J location 




ca I<,m-1 
ts lom 


^ moves 1 f^ « (m«l)^ into i(t«-mT) 
J location and loses T(t-mT) 





e 
o 






ca loO 

ts 1,1 


^ moves i(t) into i(t ■=■ T) location 




si 

rd loO 


selects the relevant input device 
(as specified by the address section) 
and makes computer wait until device 
receives a new input value 

reads content of input device into 
i(t) location 


P«,(6m*6n+6) 


sp P,l 


returns control to beginning of 

whole program, | 
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The manipulations are seen to cost 2 instructions per term of the regression 
equation * 

There are various ways in which this program can be streamlined* 
The main considerations are storage and time. It is piossible to save 
substantial storage (with a sufficiently long regression equation) by 
p370gramming the 6 instructions (k arithmetic and 2 manipulative) required 
for each term only once and using them over and over for the various terms, 
each timeo In order to do so, a short program must be added to change the 
appropriate address sections in the 6 instructions, which can thus be 
made to compute a different term each time. This address-changing routine 
materially lengthens the time of calculation, unless some very specialized 
instructions or equipment is designed. 

It appears more desirable to concentrate on reducing the time 
requirements in most control applications., for storage is easier to increase 
than speed, irfiich seems to be the ultimate limitation in the applicability 
of digital computers to controlling. In our present example a notable 
reduction in time, and also in storage, results from mixing the arithmetic 
and manipulative steps and using a new instruction, ex, which exchanges 
the content of the storage register specified by its address with the 
content of the accumiilator. The corresponding coded program, which still 
uses the same constant and data storage, follows. 



This instruction is actually used in Whirlwind applications on a temporary 
basis* The code used for this instruction is £e to indicate its temporary 
naturei final adoption of this instruction, however, is likely* 
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SECOND 


PROGRAM, 


COMPOSITE 


Register 

Noo 


Content 
(Instructions) 


Results 


P,l 


ca 0*n 








P.2 


mr B.n 








P.3 


!5X 0*n-l 


/ 


-> 


into Storage: -b %{X - nT), 
partial result " 
into AC: ^(t - (n-l)T] 

puts 'o t - (n-l)T] into o(t - nT) 
location for nex^ san?>ling.. 
AC still holds -8r(t - (n-l)f) 


P,Ii 


ts O.n 


< 


V 






P»5 


mr B«n-1 








P«6 


ad O.n-1 








Pe7 


ex 0.nT2 




^w 


into Storage: partial result 
into AC: S^lt - {n-2)T] 

•?j[t - (n-2)T] toTjf - Cn-l)Tj 

location 


y 




P*8 

* 


ts O.n-1 

• 


\ 


V 




^ 


• 
P.Hn-7 


• 
mr B.2 








etc. 


ad 0.2 










ex 0.1 


y' 


-> 


into Storage: partial result 
into AC: '^ft - T) 




ts 0.2 




-> 


ott - T) to ?yCt - 2T) location 




mr B.l 










ad 0.1 










ex I.m 


^ 


> 


into Storagei partial result 

> b -^rCt - kT) 
k=l ^ 

into AC: T(t - mT) 




mr A*m 










ad r.m 










ex I.ra-1 


*r 


■> 


into Storage: partial result 
into AC: i [t - (m-l)t| 


P*iin+3 


ts I. in 




-> 


r|t - (m-l)5 toTCt-BiT) location 
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Register 
' No. 

• 


Content 
(Instructions) 

• 


Result 


• 
P.(iin*l+nj-8) 


• 
♦ 

mr A,2 

ad 1.2 






ex I.l — 


^ into Storage: partial result 
into AC: ift - T) 




< 


-^ TCt - T) into itt - 2T) location 








rar A.l 




* 


ad 1.1 






ts 0»1 


.^ in+<^ St.f>-pflg<»2 ps^rtif^i rpsiiltr 




(note content of 0.1 has already 
been used so that this register 
is available) 




si 


selects the relevant input device 
(as specified by the address section) 
and makes computer wait until device 
receives a new input value 




rd 1.1 


reads content of input device, 
i(t) into itt - T) location 




mr A.O 






ad 0.1 






ts 0.1 


_^ -fntn StnrflgftS final r'=*sn"l t 




I^Ct) into Z{% - T) location 




si 


selects the relevant output device 
(as specified by the address section) 




re 0.1 


records output, "ott), into output 
device 


P,(Im+W6) 


sp P.l 


returns control t;o beginning of 
program 
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The above composite program is seen to result in considerable 
saving of storage and time over the first program, -which was given mainly 
for illustrative purposes. It uses four instructions per term calculated 
rather than 6, and even saves two data registers, I»0 and R»0» Register 1*0 
is not needed because the incoming data is immediately used in the calculation 
while register R,0 is superfluous because the partial results can be stored 
in the register from which the data has just been removed for calculation* 

One should note another important advantage of the second program: 
to all practical extent, it eliminates computational delays entirely. 
This is so, because all the computation is performed in advance of receiving 
the input, and when the input value i(t) is received, there are only a few 
instructions to be carried out in order to obtain the output, '^'(t)* Only 
direct regression programming has this advantage. 

The tally of direct regression composite programming in terms of 

m and n, the degree of numerator and denominator polynomials of the program 

transfer function, is as follows: 

Time requirement (in number of 
instructions to be carried out 
in sequence at each sampling) Urn * kn * 6 



Storage Requirements: 

Constants m + n ♦ 1 

Data m ♦ n 

Instruction . . i . ^ 

Total lm*hn*6 

6m ♦ 6n » 7 



The above tally is made under the assumption that none of the constants 
are zero. If some constants are zero, the constant and program storage, 
as well as the time, requirements will be reduced, but not the data storage 
requirement. These more specific requirements are taken into account in 
the summary of Art. 3.35« 
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APPENDIX C 
Coding of Cascaded Programs 

The set of regression equations 
o-^Ct) « i(t) -d^o'^(t-T) 

©^(t) = o'^(t) ^ C20^(t.T) -d202(t-T) 



(C-1) 



o(t) « 8.V^ At) * c 5-^ -(t-T)^m-T)] 
o|_ n-x n n-i n J 



is to be programmed* Assume the following arrangements of number storage: 



Register 
No. 


Content 
(Constants) 


Register 
No. 


Content 
(Data) 


D.l 
D.2 

« 

• 

D«n 


• 

« 

-d 
n 


0.1 
0.2 

. 

O.n 


^^(t-T) 
^'^(t-T) 

. 

• 
o(t - T) 


Co2 

• 
« 

C*n 
A.0 


°2 

c 

n 

a 




R.O 


Partial Result 
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In the coded program to follow it is assumed that none of the indicated 
c is zero| i»e», ra « n - 1» Variations are easily accoiinted for. The 
program inst motions follow. 



Register 
No. 


Content 
(Instruction) 


Description 


P.l 

P. 2 

P.3 
P.li 
P.5 


si 

rd R,0 
ca 0.1 
rar D.l 
ad R.O 


selects input device and waits until 
device has new input valu% ±(t) 

_^ reads i(t) into temporary location 

-d.^: (t - T) obtained 
w(t) obtained 


P.6 

F.7 
P.8 

P.9 
P.IO 
P.ll 
P«12 


ex 0.1 - 

mr C.2 

ad 0.1 

ts R.O — 

ca 0.2 

mr D.2 

ad R.O 


-> to Storage: 'SlCt) into'o'^(t-T) 
location 

— to AC: Oj^tt - T) 

Oi^tt) ♦ CgCit - T) obtained 
_^ to Storage: partial result 

'd'^it) obtained 


• 
• 

a 


• * 
• 
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Contimieds 



Register 
No. 


Content 
(Instruction) 


Description 


P.(7n-S) 
etc* 


ex O.n-1 — 

mr O.n 

ad O.n-1 

ts R.O _ 

ca 0*n 

rar D.n 

ad R.O 

mr A.O 

ts O.n - 


-^ to Storage: TS ^(t) into 

■3^ .(t-T) location 
n«-l 

— to AC: Z ,(t-T) 
n-J. 

•S - (t) ♦ c ^ T (t-T) obtained 
n-1 n n-1 

— ^ to Storage: partial result 

©"(t) obtained 

-_^ to Storage: "©(t) into '^( t-T) 
location 




si 

re O.n 


select output device 

records output, ^(t), into output 
device 


P.(7nt3) 


sp P.l 


returns control to beginning of 
program 



Thus^ if m = n - 1, the program is 7n ♦ 2 instructions long. Suppose 
m = n - 2 and let c« = 0^ then the sequence ^,6 through P.12 above would 
be replaced by the following shorter sequence P' .9 through P' .12 . 




puts "S^Ct) irito "^^(t-T) location 



'op(t) obtained 
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Thus, each c, « saves 3 instructions. 

The tally for cascade programming can now be written: 
Time Requirements? 3ni ♦ Itn ♦ ^ 

Storage Requirements 

Constants m ♦ n ♦ 1 

Data n 

Temporary 1 

Instruction 3m * Un ♦ 6 

Total Urn ♦ 6n ♦ 8 



Comparison of these requirements with those of other methods of programming 
is done in Art. 3«3^» 
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APPENDIX D 
Coding of Parallel Programs 

The set of regression equations 
Oj^Ct) - f J(t) . d^o^Ct - T) 

-o^M . fgTtt) . d^-r^Ct - T) 



(D-l) 



nsjit) . f tCt) . d T (t - T) 

n n n n 
o(t) « o-^Ct) ♦'o^W ♦ ..• ♦^(t) 

is to be programmed. Assume the following arrangement of number storages 



Register 
No. 


Content 
(Constants) 


Register 
No. 


Content 
Data 


D.l 


-^ 


6.1 


o^(t - T) 


' D.2 


-^2 


0.2 


o^gCt - T) 


• 


• 


• 


• 


• 


• 


• 


. 


Don 


-^n 


O.n 


r (t - T) 


F.l 


% 


I.O 


i(t); also o(t) 


J\2 


'2 









• 


* 




•- 


•' 






F,n 


f 

n 







None of the constants can be zero* The program instructions follow. 
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Segister 
No. 



Content 
(Instructions) 



Description 



Pa 



P.2 



si 



rd I.O 



selects input device and waits 
until device has new input value; 

T(t) 

reads t(t) into its assigned 
storage register 



P.3 
P.I4 
P.^ 

P.6 
P.7 
P.8 



ca 1*0 
mr F.l 
ex 0.1 

mr D,l 
ad 0.1 
ts 0.1 



to Storage J f JiTt) 
to AC J ^^(t - T) 



to Storage: 3", (t) 



P.9 

P.IO 

P.ll 

P.12 
P.13 
P.lU 



ca 1.0 
mr F.2 
ex 0.2 

mr D.2 
ad 0.2 
ts 0.2 



to Storage: f^itt) 
to AC: ^^(t - T) 



to Storage: ol(t) 



P.C6n-3) 
etc. 



ca 1.0 
mr F.n 
ex O.n 



■> to Storage: f T(t) 

r " n 

- to AC: cT Ct - T) 
n 
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Register 
No* 


Content 
(Instructions) 


Description 




rar D.n 
ad O.n 

+.«? O.n 


X- +.rt S+.nT*nor<ai r\ ^ + ^ 






P.(6n+3) 
etc. 

e 
a 
» 

P.(7n+2) 


ad 0*n-l 

ad 0,nw2 

• 
• 
• 

ad 0.2 
ad 0.1 
ts 1*0 


o'Ct)*o' ,(t) obtained 
n n-1 

etc. 

'o(t) obtained 
-^ to Storages lS(t) 


etc. 


si 

re I.O 

sp P.l 


selects output device 

records output, "oCt) into 
output device 

returns control to beginning of 

program 



In parallel programming none of the indicated constants can be 
seroj and the only possible saving is when several constants have the same 
value. Even then the program itself is not affected materially. 

The tally for parallel programming follows s 

Time Requirement yn ♦ 5 

Storage Requirements s 

Constants 

Data 

Instructions 

Total 





2n 






n 


4' 


1 


7n 


+ 


^ 


lOn 


♦ 


6 





Further discussion of these requirements is left to Art. 3«35» 



